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^ (57) Abstract: Methods, apparatus and computer program products provide efficient techniques for reconstructing surfaces from 
^ data point sets. These techniques include reconstructing surfaces from sets of scanned data points that have preferably undergone 
— . preprocessing operations to improve their quality by, for example, reducing noise and removing outliers. These techniques include 
£5 reconstructing a dense and locally two-dimensionally distributed 3D point set (e.g., point cloud) by merging stars in two-dimensional 
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Q plurality of points pi in a 3D point set S that at least partially describes the #D surface, by projecting the plurality of point pi onto 

planes Ti that arc each estimated to be tangent about a respective one of the plurality of points pi. The plurality of stars arc then 
^ merged into a digital model of the 3D surface. 
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METHODS, APPARATUS AND COMPUTER PROGRAM PRODUCTS 
THAT RECONSTRUCT SURFACES FROM DATA POINT SETS 

Reference to Priority Application 
This application claims priority to U.S. Application Serial No. 
10/152,444, filed March 21, 2002 and U.S. Provisional Application Serial 
No. 60/324,403, filed September 24, 2001 , the disclosures of which are 
5 hereby incorporated herein by reference. 

Field of the Invention 
This invention relates to methods and systems that reconstruct 
three-dimensional (3D) surfaces and, more particularly, to methods and 
systems that reconstruct 3D surfaces from data point sets. 
10 Background of the Invention 

Conventional techniques that use differential concepts to reconstruct 
surfaces from scanned data point sets typically assume a finite set of 
points that are sampled on the surface of a shape in a three-dimensional 
space and ask for an approximation of that surface. Such techniques may 
15 be classified by the assumptions they make about the data point sets. 

Techniques that make structural assumptions may simplify the 
reconstruction task by providing the points in a specific order. Techniques 
that make density assumptions may enable the application of differential 
concepts by providing sufficiently many point samples within each data 
20 point set. Techniques that avoid assumptions typically require the 

reconstruction operations to rely on general principles of describing 
geometric shapes. 
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One conventional technique that incorporates density assumptions 
is described in an article by H. Hoppe et at., entitled "Surface 
Reconstruction from Unorganized Points," Computer Graphics, 
Proceedings of SIGGRAPH, pp. 71-78 (1992). This technique uses normal 
5 . estimates to generate a signed distance function from which a surface is 
. extracted as a zero-set Another conventional technique is described in an 
article by N. Amenta et at., entitled "Surface Reconstruction by Voronoi 
Filtering," Discrete Computer Geometry, Vol. 22, pp. 481-504 (1999). This 
technique exploits shape properties of three-dimensional Voronoi cells for 
. 10 densely sampled data points. Unfortunately, because these reconstruction 
techniques rely heavily on the quality of the data, they may fail if a given 
data point set does not adequately support the application of differential 
concepts. For example, these reconstruction techniques may fail if there 
are large gaps in the distribution of the data points in a set or if the 
15 accuracy of the data points is compromised by random noise. These 

techniques may also fail if the data point sets are contaminated by outliers. 
Data point sets having relatively large gaps, high levels of random noise 
and/or outliers may result from scanning objects having sharp edges and 
comers. 

20 Additional surface reconstruction techniques can be distinguished 

based on the internal operations they perform. For example, in "sliced 
data" reconstruction techniques, the data points and their ordering are 
assumed to identify polygonal cross-sections in a finite sequence of parallel 
planes. This assumption may simplify the complexity of the technique, but 

25 it also typically limits the technique to data generated by a subclass of 

scanners. A survey of work that includes this technique is described in an 
article by D. Meyers et al., entitled "Surfaces from Contours," ACM Trans, 
on Graphics, Vol. 1 1, pp. 228-258 (1992). In another feconstruction 
technique, the data points are used to construct a map f:R z - R, and the 

30 surface is constructed as the zero set, V 1 (0). The zero set may be 

constructed using a marching cube algorithm that is described in an article ' 
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by W. Lorensen et al., entitled "Marching Cubes: A High Resolution 3D 
Surface Construction Algorithm," Computer Graphics. Proceedings of 
SIGGRAPH, Vol. 21, pp. 163-169 (1987). One example of this 
reconstruction technique is described in the aforementioned article by H. 
5 Hoppeetal. Attempts to generalize and apply surface meshing techniques 
to the reconstruction of surfaces from unstructured data point sets have 
also been presented in a survey paper by S.K. Lodha and R. Franke, 
entitled -Scattered Data Techniques for Surfaces," Proceedings of a ' 
Dagstuhl Seminar, Scientific Visualization Dagstuhl '97, Hagen, Nieison 
and Post(eds.), pp. 189-230. Additional techniques for automatically 
wrapping data point sets into digital models of surfaces are also disclosed 
in U.S. Patent No. 6,377.865 to Edelsbrunner et al., entitled "Methods of 
Generating Three-Dimensional Digital Models of Objects by Wrapping 
Point Cloud Data Points," assigned to the present assignee, the disclosure 
15 of which is hereby incorporated herein by reference. 

Summary of the Invention 
Methods of modeling three-dimensional (3D) surfaces include 
preferred techniques to reconstruct surfaces from respective sets of data 
points that at least partially describe the surfaces. The data points may be 
generated as point clouds by scanning an object in three dimensions. 
These reconstruction techniques may include improving the quality of the 
reconstructed surfaces using numerical approximations of local differential 
structure to reduce noise and remove outliers from the data points. The 
numerical approximations of local differential structure may be achieved 
using differential concepts that include determining surface normals at 
respective ones of the data points. A local differential structure of a bundle 
of the surface normals can be used to define principal curvatures and their 
directions. In particular, a number of substantially non-zero principal 
curvatures may be used to determine the types of approximating surfaces 
that can be locally fit to the data points. By construction, these 
approximating surfaces are typically only sensitive to an average amount of 



20 
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, - . noise that is locally present in the data points. Accordingly, the 
approximating surfaces vary considerably slower than the points 
themselves. This difference is preferably exploited in reducing the local 
variability of the data points and, therefore, the amount of random noise 
5 present therein. The same differential concepts may also be used to detect 
outliers and to subsample the data points in a curvature sensitive manner. 

Methods of modeling 3D surfaces according to first embodiments of 
the present invention include determining an estimated normal for each of a 
plurality of points in a 3D data point set that at least partially describes the 
10 3D surface. The plurality of points may constitute all or less than all of the 
points in the 3D point set A differential structure of the estimated normals 
is evaluated to estimate principal curvature directions on the 3D surface 
and to classify a respective local neighborhood of each of the plurality of 
points in terms of its shape characteristic. An approximating surface is 
15 then determined for each of the local neighborhoods. A denoising 

operation may be performed on the 3D point set by moving each of the 
plurality of points to a respective approximating surface that is associated 
with a local neighborhood of the respective point. Other methods of 
modeling 3D surfaces may include operations to determine a respective set 
20 of near points S, for each of a plurality of points p, in a 3D point set S that at 
least partially describes the 3D surface, where 5; c S. An operation may 
then be performed to determine a normal bundle for the 3D point set S by 
determining a respective plane ft, of best fit for each of the sets of near 
points S, and a respective normal r>, for each of the planes ft, of best fit. A 
25 respective approximating surface is then determined for each of the sets of 
near points S { using the normal bundle to estimate respective principal 
curvature directions for each of the sets of near points S ( . 

Methods of modeling 3D surfaces according to still further 
embodiments of the present invention include determining a respective set 
30 of near points for each of a plurality of points in a 3D point set that at least 
partially describes the 3D surface and then fitting each set of near points 
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with a respective approximating surface. The approximating surfaces may 
be selected from the group consisting of planes, cylinders and quadratic or 
cubic surfaces. The 3D point set is then denoised by moving each of the 
plurality of points in the 3D point set to the approximating surface 
5 associated with its respective set of near points. Additional embodiments 
include denoising a 3D point set that at least partially describes the 3D 
surface by classifying a first neighborhood of points in the 3D point set 
using (i) a mass distribution matrix (MDM) of the first neighborhood of 
points to estimate first normals associated with the first neighborhood of 

1 0 points and (ii) a normal distribution matrix (NDM) of the first normals to 

estimate principal curvature directions. Operations are then performed to 
fit an approximating surface to the first neighborhood of points and then 
move at least one point in the first neighborhood of points to the 
approximating surface to thereby reduce noise in the first neighborhood of 

15 points. 

Additional methods of modeling three-dimensional (3D) surfaces 
include techniques to reconstruct surfaces from sets of scanned data 
points that have preferably undergone preprocessing operations to improve 
their quality by, for example, reducing noise and removing outliers as 

20 described above. These, methods preferably include reconstructing a - 
dense and locally two-dimensionally distributed 3D point set (e.g., point 
cloud) by merging stars in two-dimensional Delaunay triangulations within 
estimated tangent planes. These two-dimensional Delaunay triangulations 
may have vertices with non-zero weights and; therefore, may be described 

J5 as weighted Delaunay triangulations. In particular, these methods include 
determining a plurality of stars from a plurality of points p, in a 3D point set 
S that at least partially describes the 3D surface, by projecting the plurality 
of points p, onto planes T, that are each estimated to be tangent about a 
respective one of the plurality of points p t The plurality of stars are then 

•0 merged into a digital model of the 3D surface. 
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. . , - . The operations to determine a plurality of stars preferably include . 
. . identifying a respective subset of near points .S,- for each of the plurality of 
points Pi and projecting a plurality of points p $ in each subset of near points 
S, to a respective estimated tangent plane T t In particular, for each of a 
5 plurality of. estimated tangent planes, 7„ a star of the projection of a 

respective point p, onto the estimated tangent plane 7] is determined. The 
star of the projection of a respective point p, onto the estimated tangent 
* plane T, constitutes a two-dimensional (2D) Delaunay triangulation (e.g., 
2D weighted. Delaunay triangulation). To improve efficiency, the operation 
10 to identify a respective subset of near points S, for each of the plurality of 
: points & may include storing the 3D point set S jn an oct-tree. This 
operation may also include determining a width 2r 0 of a near point search 
cube using a random sample RcS, where r 0 is a positive real number. 
Then, for each of the plurality of points p„ a subset of ^ points that are 
15 closest in Euclidean distance to p x may be determined, where is a 

positive integer. From this subset, all closest points that are also within a 
respective near point search cube, which is centered about a 
corresponding point p, and has a width equal to 2r 0 , are selected. 

Additional techniques to reconstruct surfaces include modeling a 
20 threenJimensional (3D) surface by determining a plurality of stars from a 
plurality of points p t in a 3D point set S that at least partially describes the 
3D surface, by projecting each of the plurality of points p t onto a respective 
plane TJ that is estimated to be tangent about the corresponding point p x . 
Weights, which are based on projection distance (e.g., (projection 
25 distance) 2 ), are also preferably assigned to each projected point so that 
subsequent operations to merge triangles and edges result in fewer 
conflicts and, therefore, fewer holes in resulting digital models. The 
plurality of stars are then merged into a model of the 3D surface. The 
merging operation includes eliminating edges and triangles from the 
30 plurality of stars that are in conflict and merging nonconflicting edges and 
triangles into a surface triangulation. The merging operation may include 
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sorting triangles within the plurality of stars and removing those sorted 
triangles that are not in triplicate. The sorted triangles that have not be 
removed are then connected to define a triangulated pseudomanifold as a 
two-dimensional simplicial complex in which edges and triangles of a star 
that share a vertex form a portion of an open disk. Operations may also be 
performed to sort edges within the plurality of stars that do not belong to 
any of the triangles in the triangulated pseudomanifold and remove those 
sorted edges that are not in duplicate. The sorted edges that have not 
been removed may then be added to the triangulated pseudomanifold. 

Still further methods include improving the quality of a scanned data 
point set by removing noise and/or eliminating outliers therefrom and/or 
subsampling the point set and then reconstructing a surface from the 
improved data point set. These methods preferably include determining a 
respective set of near points S t for each of a first plurality of points pf, in a 
first 3D point set S1 that at least partially describes the 3D surface, where 
S, c S1. Each set of near points § is fit with a respective approximating 
surface. A denoising operation is then performed which converts the first 
3D point set S1 into a second 3D point set S2. This denoising operation is 
performed by moving at least some of the first plurality of points pf, in the 
first 3D point set S1 to the approximating surfaces associated with their 
respective sets of near points S,. Although rare, If the operation to move a 
point to a respective approximating surface requires a spatial translation 
that exceeds a designated threshold value, then the point may be removed 
as an outlier. A plurality of stars are then determined from a second 
plurality of points p2, in the second 3D point set S2. These stars are 
determined by projecting the second plurality of points p2, onto planes 7; 
that are estimated to be tangent about a respective one of the second 
plurality of points p2 h The stars are then merged into a digital model of the 
3D surface. 

The operations to fit each set of near points S, with a respective 
approximating surface comprise fitting a first set of near points S 1 with a 
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first approximating surface by determining respective planes /ij of best fit for 
each of a plurality of points % in the first set of near points S, and then 
determining an estimated normal n $ for each of the points ^ as a. normal of 
its respective plane h { of best fit. A shape characteristic of the first set of 
5 near points S 1 is also classified by determining estimates of principal 

curvature directions for a point pf, from a plurality of the estimated normals 
n Y The shape characteristic may be plane-like and/or edge-like and/or 
comer-like. 

According to still further preferred aspects of these embodiments, 

10 the operations to determine a plurality of stars include operations to 
determine weighted Delaunay triangulations on tangent planes. In 
particular, these operations may include determining a star of a, first point in 
a 3D point set that at least partially describes the 3D surface, by: (i) 
projecting the first point and at least second, third and fourth points in a 

1 5 neighborhood of the first point in the 3D point set to a plane, (ii) assigning . 
respective weights to each of the second, third and fourth points that are 
based on projection distance, and (iii) evaluating whether the projection of 
the fourth point is closer than orthogonal to an orthocenter of a first triangle 
having projections of the first second and third points as vertices. If the 

20 projection of the fourth point is closer than orthogonal, then It Is included as 
a vertex of a new triangle in the star of the projection of the first point and 
the first triangle is discarded. If the projection of the fourth point is farther 
than orthogonal, then it is included as a vertex of a new triangle containing 
the first, third and fourth points as vertices. 

25 Additional embodiments of the present invention may also include 

operations to create stars of projected points by sequentially connecting a 
neighborhood of projected points on a plane to a first projected point on the 
plane by evaluating whether at least one projected point in the 
neighborhood of projected points is closer than orthogonal to an 

30 orthocenter of an orthocircle. This orthocircle is defined by a triangle 
containing the first projected point and two projected points in the 
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. neighborhood of projected points as vertices, with the neighborhood of 
* projected points having weights associated therewith. These weights are 
each a function of a projection distance between a respective projected 
point in the neighborhood of projected points and a corresponding point in 
a 3D point set that at least partially describes the 3D surface. 

Still further operations may include projecting a first point in a 3D 
point set that at least partially describes the surface and a set of points in a 
neighborhood of the first point to a plane that is estimated to be tangent to 
the surface at the first point An operation may then be performed to create 
a weighted Delaunay triangulation comprising triangles that share a 
projection of the first point in the plane as a vertex and include at least 
some of the projections of the set of points in the neighborhood of the first 
point as vertices that are weighted as a function projection distance 
squared. 

Embodiments of the present invention may also include operations 
to model a three-dimensional (3D) surface by projecting a first point in a 3D 
point set that at least partially describes the surface and a set of po(nts in a 
neighborhood of the first point to a plane that is estimated to be tangent to 
the surface at the first point Additional operations include creating a 
weighted Delaunay triangulation comprising triangles that share a 
projection of the first point in the plane as a vertex and include at least 
some of the projections of the set of points in the neighborhood of the first 
point as vertices that are weighted as a function of projection distance. 
Operations to create a weighted Delaunay triangulation may include 
operations to evaluate whether or not one or more of the projections of the 
set of points in the neighborhood of the first point are closer than 
orthogonal to an orthocenter of a first triangle in the weighted Delaunay 
triangulation. The operations to create a weighted Delaunay triangulation 
may include evaluating a 4x4 matrix containing coordinates of the vertices 
of the first triangle as entries therein along with entries that are functionally 
dependent on the weights associated with the vertices of the first triangle. 
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v , * These operations to evaluate the matrix may include computing a 
determinant of the matrix 

Brief Description of the Drawings 
FIG. 1 illustrates a closed axis-aligned cube that may be used when 
5 determining collections of near points. 

FIG. 2 illustrates a standard triangle A 1 A 2 A 3 -wtth shading that 
corresponds to data point sets that are strongly comer-like, edge-like and 
plane-like. 

FIG. 3 illustrates construction of an ocl-tree in two dimensions. 
10 FIG. 4 illustrates shaded stars within a weighted Delaunay 

triangulation. 

FIG. 5 Illustrates sequential operations to construct a star. 

FIG. 6 illustrates a trist data structure comprising an array of 
(unordered) triangles connected to each other by adjacency. 
1 5 FIG. 7 illustrates five dummy triangles, four of which connect co to 

boundary edges and one of which connects co to a principal edge in a star 
ofp/. 

FIG. 8 illustrates two 'crossing' principal edges and two 'crossing' 
triangles. 

20 FIG. 9 illustrates a hole whose cycle of seven directed edges 

includes two along a principal edge of a pseudomanifold, 

FIG. 1 0 illustrates an embedding of a directed boundary cycle that 
cuts a plane into an open disk on its left and one or more shaded open 
regions on its rights. 

25 FIG. 1 1 illustrates a set of points, including p 1 and p 2 , that define a 

portion of a three-dimensional (3D) surface. 
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FIGS. 1 2A-12D are illustrations that describe operations to create 
stars by sequentially connecting points projected to respective tangent 
planes and generate stars, according to embodiments of the present 
invention. 

5 FIG. 13 illustrates a 3D triangulation that results from merging the 

stars illustrated by FIG. 12D. 

FIGS. 14A-14G illustrate operations to sequentially connect . 
projected points according to embodiments of the present invention. 
FIG. 1 5 is a flow-diagram that illustrates operations for modeling 
1 0 three-dimensional surfaces according to embodiments of the present 

invention. 

FIG. 16 is a flow-diagram that illustrates operations for modeling 
three-dimensional surfaces according to embodiments of the present 
invention. 

15 FIG. 17 is a flow-diagram that illustrates preferred operations for 

merging stars according to embodiments of the present invention. 

FIG. 18 is a flow-diagram that illustrates operations for 
preprocessing three-dimensional point sets to reduce noise, according to 
embodiments of the present invention. 
20 FIG. 1 9 is a flow-diagram that Illustrates operations for 

preprocessing three-dimensional point sets to reduce noise, according to 
embodiments of the present invention. 

FIG. 20 is a flow-diagram that illustrates operations for modeling 
three-dimensional surfaces according to embodiments of the present 
25 invention. . 

FIG. 21 is a block diagram that illustrates a general hardware 
description of a computer workstation that performs operations according . 
to embodiments of the present invention. 
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3*£<ty. ' Description of a Preferred Embodiment 1 

The present invention now will be described more fully hereinafter 
with reference to the accompanying drawings, in which preferred 
embodiments of the invention are shown. This invention may, however, be 
5 embodied in many different forms and applied to other articles and should 
not be construed as limited to the embodiments set forth herein; rather, 
these embodiments are provided so that this disclosure will be thorough 
and complete, and will fully convey the scope of the invention to those 
. skilled in the art. The operations of the present invention, as described 
10 more fully hereinbelow and in the accompanying figures, may be performed 
by an entirely hardware embodiment or, more preferably, by an 
embodiment combining both software and hardware aspects and some 
degree of user input. Furthermore, aspects of the present invention may 
take the form of a computer program product on a computer-readable 
15 storage medium having computer-readable program code embodied in the 
medium. Any suitable computer-readable medium may be utilized 
including hard disks, CD-ROMs or other optical or magnetic storage 
devices. Like numbers refer to like elements throughout 

Various aspects of the present invention are illustrated in detail in 
20 the following figures, including flowchart illustrations. It will be understood 
that each of a plurality of blocks of the flowchart illustrations, and 
combinations of blocks in the flowchart illustrations, can be implemented by 
computer program instructions. These computer program instructions may 
be provided to a processor or other programmable data processing 
25 apparatus to produce a machine, such that the instructions which execute 
on the processor or other programmable data processing apparatus create 
means for implementing the operations specified in the flowchart block or 
blocks. These computer program instructions may also be stored in a 
computer-readable memory that can direct a processor or other 
30 programmable data processing apparatus to function in a particular 

manner, such that the instructions stored in the computer-readable memory 
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produce an article of manufacture including instruction means which 
implement the operations specified in the flowchart block or blocks. 

Accordingly, blocks of the flowchart illustrations support 
combinations of means for performing the specified operations, 
combinations of steps for performing the specified operations and program 
instruction means for performing the specified operations. It will also be 
understood that each of a plurality of blocks of the flowchart illustrations, 
and combinations of blocks in the flowchart illustrations, can be 
implemented by special purpose hardware-based computer systems that 
perform the specified operations or steps, or by combinations of special 
purpose hardware and computer instructions. 

Methods of modeling three-dimensional (3D) surfaces according to 
embodiments of the present invention include preferred techniques to 
reconstruct a surface from data points that at least partially describe the 
surface. The data points may be generated as point clouds by scanning an 
object in three dimensions. These reconstruction techniques may include 
operations to improve the quality of the reconstructed surfaces using 
numerical approximations of local differential structure to reduce noise and 
remove outliers from trie data points. The numerical approximations of 
local differential structure may be achieved using differential concepts that 
include determining surface normals at respective ones of the data points. 
A local differential structure of a bundle of the surface normals can be used 
to define principal curvatures and their directions. In particular, a number 
of substantially non-zero principal curvatures may be used to determine the 
types of approximating surfaces that can be iocally fit to the data points. 
By construction, these approximating surfaces are typically only sensitive to 
an average amount of noise that is locally present in the data points. 
Accordingly, the approximating surfaces vary considerably slower than the 
points themselves. This difference is preferably exploited in reducing the 
local variability of the data points and, therefore, the amount of random 
noise present therein. The same differential concepts may also be used to 
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; • detect outliers and to subsample the data points in a curvature sensitive 
manner. 

The numerical approximations of local differential structure are 
preferably obtained using mean square error minimization techniques. A 
5 mass distribution matrix (MDM) of a local collection of data points is used 
to estimate normals and a normal distribution matrix (IMDM) of a local 
collection of normals is used to estimate principal curvature directions. In 
both cases, a preferred numerical method includes the spectral 
decomposition of the matrix, which provides direction estimates through 
10 eigenvectors and variability estimates through eigenvalues. Numerical 
stability can be improved by normalizing all input data and by using 
standard operations for the singular value decomposition of the matrices. 
Because a time consuming operation includes finding a set of near points 
for each of a plurality of points in the data point set, an implict oct-tree is 

15 used to store the data point set in an efficient manner. 

In particular, operations to improve the quality of surfaces generated 
by surface reconstruction techniques include identifying collections S,for 
each of a plurality points p, (and possibly all points) in a finite three- 
dimensional (30) set S of points that are relatively densely distributed on a 

20 surface in R 3 , where S, c S and p, e Each collection S, is referred to 
herein as the set of "near points" of p, and is used as a primary source of 
information about the neighborhood of p, on the surface. In particular, S, is 
the collection of points p y eS that He within a closed and axis^aHgned cube 
10 that has a side-length of 2r 0 and is centered at p t This is illustrated by 

25 FIG. 1, where the solid black points located within the cube 10 belong to S t 
The real number constant r 0 is chosen such that an expected size of the 
sets of near points S f is some constant yield, which may be a user-defined 
integer. Typical values for low, medium and high yield are 25, 50 and 100, 
respectively. As described more fully hereinbelow, to compute a set of near 

30 points Sf for each of a plurality of points p h the data point set S c R 3 is stored 
in an implicit oct-tree. A random sample R c S is used to compute a width 
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2r 0 of a closed and axis-aligned cube 10 having a specified (e.g., user 
specified) average yield. Then, for each point p, e S, a subset of ^ nearest 
points is computed. Those points in the subset that also lie within the cube 
of width 2r„ (centered at p,) are then selected as points in a respective near 
5 point set S,. The validity of the operations described herein are independent 
of the value of the specified average yield, but the value of the yield 
influences the results they generate and the computing time required to 
generate the results. 

For a given set of near points S, corresponding to a point p, e S, the 

1 0 plane that minimizes a mean square distance to the near points necessarily 
passes through the centroid of the points. The near points may be 
translated by minus the centroid or, equivalents, the centroid can be treated 
as the origin and only planes of zero offset relative to the origin need be 
considered. Each such plane can be written as the set of points xeTSL 3 that 

1 5 satisfy the relationship x T u=0 for some unit normal u e S 2 . The sum of the 
square distances of the points p 7 e S,from a plane having zero offset is: 

£»=2>J«) 2 

J 
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The mass distribution matrix M, <* YPPi is symmetric and positive semi- 
definite and thus has three non-negative real eigenvalues //, * ^ 
Accordingly, M, = ULU T , where: 



25 
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The columns of U are the corresponding unit eigenvectors e lt ej and e 3 in 
the same order. Multiplication with If puts a point p y into the coordinate 
frame spanned by the eigenvectors, and multiplication with t/puts it back 
into the original coordinate frame. For a direction u = « l e 1 + + w 3 e 3 , the 
sum of square distance is E£u) = //,("i 2 ) + P&z) + th( u *Y The preimage of 
unity is the ellipsoid Ef 1 <1). and the half-axes of £," 1 (1) have lengths equal to 
(//,)"* along u { . The unit vector with smallest error is therefore in th£ direction 
of the longest half-axis, which is parallel to e 3 . The plane of best fit, h f , 
therefore consists of the points x that satisfy x T e 3 = 0. The plane of best fit, 
h i9 typically does not pass through p h The unit normal /?, to the plane of best 
fit ^ is treated herein as the estimated normal at p h To improve run time for 
situations involving three-by-three matrices, which are the most common, 
the eigenvalues are computed as the roots of the characteristic polynomial. 
The first root is approximated by Newton iteration, and the other two by 
solving the remaining quadratic equation analytically. Other techniques for 
computing eigenvalues may also be used. 

The estimated normals n, at the points p, are used to estimate 
principal curvature directions on the surface (i.e., two principal curvature 
directions for each point p t ). These estimates of principal curvature 
directions are then used to classify the local neighborhood of each point P; in 
terms of its shape characteristics. In particular, the function F; R 3 - a is 
used that sums the square distances to the normal planes passing through 
the origin, which is defined by: 

^(x)=2(x r /i,) 2 
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Similar to the mass distribution matrix M„ the normal distribution matrix N, = 
= In/)/ is symmetric and positive semi-definite with non-negative 
eigenvalues v, * v 2 * v 3 . A large eigenvalue corresponds to a direction in 
which the sum of squared scalar products or. equivalents, the sum of 
square distances to the normal planes increases quickly. Accordingly, a set 
of near points S, can be treated as having plane-like configuration if there is 
only one such principal curvature direction or, equivalents, only one large 
eigenvalue. 

Operations to classify a local neighborhood of each point p, in terms 
of its shape characteristics continue by normalizing the eigenvalues to f, = 
(VxXv, + v 2 + v 3 )-\ for / = 1 , 2 and 3, where t, + 1 2 + t 3 = 1 and t, * f 2 * f 3 . As 
illustrated by FIG. 2, the triplets of numbers t, that satisfy these conditions 
define a triangle in a barycentric subdivision of a standard triangle. The 
smaller two eigenvalues correspond to the two principal curvatures. 
Neighborhoods with two, one and zero substantially norvzero principal 
curvatures are classified using two constants e 0 , e, that meet the following 
relationship: 0 < e„ < e, < v&. The shape classification thresholds e„ and e, 
may be set at 0.0075 and 0.015, respectively. A set of near points S, is 
treated as strongly plane-like, edge-like or comer-like based on the following 
relationships: ' . 

plane-like - e 0 * t 2 a f 3 
edge-like - c 0 * t s and 6,54 
comer-like - e, * t 3 s 

1 

Similarly, a set of near points S, is treated as being weakly plane-like, 
weakly edge-like and weakly comer-like if the normalized eigenvalues 
satisfy the same inequalities with e 0 and e, interchanged. As illustrated by 
FIG. 2, the points + + fyAj populate the lower left triangle in the 
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barycentric subdivision of the standard triangle 20. The light, 
medium and dark shaded areas correspond to data point sets, that are 
strongly comer-like, strongly edge-like and strongly plane-like, respectively. 
The white strips between the shaded areas correspond to data point sets 
that have at least two weak classifications. FIG. 2 also illustrates that each 
set of near points S, has at most one strong classifier, and if it has no strong 
classifiers then it has at least two weak classifiers. More precisely, the 
region of multiple weak classifications overlaps the region of strong 
classifications only along a boundary. The normalized eigenvalues can 
therefore be used to express each multiple weak classification as a linear 
combination of strong classifications. 

The normalized eigenvalue j, equals (f, - e 0 y(e n - e 0 ), provided that e 0 
s^e lf and this is the solution to the linear interpolation f, = (1- s,)e 0 + s& 
between e 0 and e v The normalized eigenvalues s 2 and s 3 are used to 
linearly interpolate between strong classifications. For Case 1, where f 3 < e 0 
* t 2 s the set of near points S, is both weakly plane-like and weakly edge- 
like, but not weakly comer-like. The normalized eigenvalue s 2 is defined, 
and for S, the fraction 1-5 2 is considered strongly plane-like and the fraction 
j 2 is considered strongly edge-like. For Case 2, where e 0 ^ 3 ^ i 2 s e„ the 
set of near points S, has ail three weak classifications. For the set of near 
points S„ the fraction 1-$ 2 -s 3 is considered plane-like, the fraction s 2 is 
considered edge-like and the fraction j 3 is considered corner-like. For Case 
3, where e 0 < t 2 < < t 2t the set of near points S,- is both weakly edge-like 
and weakly comeMike, but not weakly plane-like. For the set of near points 
S„ the fraction 1-j 3 is considered strongly edge-like and the fraction s 2 is 
considered strongly comer-like. 

In regions of low curvature, the local surface can be approximated by 
a plane and the plane of best fit /?, can be determined directly from the mass 
distribution matrix (MDM). In regions of non-trivial curvature, the local 
shape can be approximated by surfaces that are more complicated than 
planes. Cylinders with conic cross-sections can be used for near point sets 
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Sf that are edge-like and paraboloids can be used for near point sets S, that , 
are comer-like. 

The set of near points S y is edge-like if the normals lie roughly along a 
great-circle of the sphere of directions. The best approximation of the 
direction normal to that great-circle is the eigenvector that corresponds to 
the smallest eigenvalue of the normal distribution matrix (NDM). This is an 
approximation of the minor principal curvature direction, which will be 
referred to herein as a fold direction. The local shape can be approximated 
by a cylinder constructed by sweeping a line parallel to the fold direction 
along a conic in the orthogonal plane. Let g, be this plane passing through 
p /t and let S/ consist of the points in S, projected along the fold direction 
onto g h The conic g, is constructed by least square optimization. The family 
of conies considered are the zero-sets of the functions G; R 2 - R defined by: 

where x 1 and xj are the coordinates in g, measured along the other two 
eigenvector directions, using p, as the origin. The family of conies contains 
pairs of intersecting lines as limits of progressively narrower hyperbolas. 
These intersecting lines are important because they model common sharp 
edges on the surface. Multiplying the entire polynomial by a constant does 
not change the zero-set, so the assumption that £ a, = 1 can be used. The 
function G, may be treated as an affine function R s - R that takes a point in 
homogeneous coordinates, x T = (x, 2 , x, x 2% x? % x 1f Xj, 1 ), to the residual, G/x). 
Each point p/ e S/ defines such a point x y e R 5 . The affine function that 
minimizes a sum of the squared residuals is: 

J S 
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where a T = (a v a 3 , a 4 , a 5l a 6 ). The six-by-six matrix X, = £xjx/ is again 
symmetric and positive semi-definite and thus has six non-negative real 
eigenvalues. The affine function of best fit is determined by the unit 
eigenvector a that corresponds to the smallest of the six eigenvalues. 

When the matrix X,, which is computed for an edge-like 
neighborhood, has two almost equally small eigenvalues, then the conic of 
best fit may be ambiguous. In this case, a preferred operation considers 
both corresponding conies and projects the point to the closest of the points 
computed for both conies. 

The sets of near points S, that are comer-like permit the largest 
amount of variation and, therefore, are typically the most difficult to locally 
approximate. Ideally, a family of surfaces that contains triplets of 
intersecting planes is used. The family of zero-sets of cubic polynomials in 
three variables contains such triplets but leads to linear systems with as 
many as twenty unknowns. Contrasting this with the fact that points with 
comer-like neighborhoods are naturally the least common in surfaces 
bounding physical artifacts, an option is pursued that allows a solution of a 
significantly smaller optimization problem. The specific problem is finding 
the optimum function in the family of functions G, defined above. The 
coordinate frame, within which the functions G, are considered, consist of 
the three unit eigenvectors of the normal distribution matrix (NDM). 
Specifically, x 3 is measured in the surface normal direction, which is 
approximated by the eigenvector that corresponds to v„ and x 2 are 
measured in the other two eigenvector directions. 

The minimization problem is different from the above because the 
graph of G„ and not the zero-set of G„ is of interest G, is again interpreted 
as an affine function S 5 - R. For each p y e S„ x, is defined as described 
above and z } is treated as the x 3 -coordinate of p, within the frame of 
eigenvectors. The new residual is G/x ; ) - z i and the sum of squared 
residuals is computed as: 
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j 

J 

where: *=2>,*y 

To compute the minimum a, all partial derivatives are set to zero. This 
5 yields: 

ZT --—JL.a + a X, — 2 v. 

= XZa+a T X a -2y a% 

for 1 * 7*6, where X,' is the /th column of X, and y a is the Ah entry in y,. Thus. 

0 the minimum a is the solution to the linear system X^a = y t . This system can 
be treated as well-conditioned. 

The purpose of computing the approximating surfaces is to improve 
the quality of the description of the surface defined by the data point set. 
The primary goals are to reduce noise and remove outliers. The type of 

5 noise that is typically reduced is manffested by points that are close to but 

extend slightly above or below the surface. It is a common phenomenon in 
scanned data and can be caused by a variety of shortcomings of the 
scanning hardware. The approach to reduce noise is based on the 
observation that the approximating surface varies only slightly for nearby 
points, simpiy because it is computed from sets of near points S,that are 
largely the same. Noise is preferably reduced by moving each point p, onto 
the approximating surface computed for Sf Three cases are distinguished 
. as follows. 

In the first case, S, is plane-like. The approximating surface is the 
plane of best fit h, with normal vector n,. By construction, the centroid of S, 
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lies on this plane. The point p, is moved to its orthogonal projection on h h 
The orthogonal projection of the point is p/Y where: 

5 and the centroid is represented as p t . The point p/' may also be referred 

to as the de-noised location of p/. 

In the second case, S,is edge-like. The surface is a cylinder whose 
cross-section in g, is the conic of best fit, G; 1 (0). The gradient of G, at a 
point xeg, is: 

10 = 

By construction; p, lies at the origin of g„ which implies that the gradient at P; 
is: 

i VGi(0) = [a 4 ,a 5 ] r , 

The point p, can be moved iteratively in small steps along the gradient. 
15 Since p, is mostly already very close to the conic, this iterative procedure is 
simplified and the point p, is moved in one. step. Formally, the projection of 
p, is computed as: 

p^t[a 49 a s ] T 9 

such that 

20 <?,(/>/) = t 2 {a x al +a 2 a 4 a s +a 3 a 3 2 ) + r(a^ +a 3 2 ) + a 6 

vanishes. Both roots of this quadratic polynomial are computed, and the 
point p," that is closest to point p ; is selected. 

In the third case, S, is comer-like. The approximating surface is the 
25 graph of the quadratic function G, of best fit. The point p, is projected in the 
x 3 -direction onto that graph. By construction, p, is the origin of the 



-22- 



03/027961 



PCT/USO2/24220 



coordinate system, so A"=( n .0.G,(0)), again expressed in the local 
coordinate frame spanned by the eigenvectors of the normal distribution 
matrix For strongly classified sets S„ the point p," is substituted for p h For 
sets S, with mixed weak classifications, two or three points p^are computed 
and the linear combination is substituted for p,. In each case, the point p," is 
dropped from the computations^ the fit of the surface is not sufficiently tight, 
or more specifically, If the normalized smallest eigenvalue exceeds some 
constant So. Such cases are unlikely for plane-like and edge-like sets, 
because the surfaces typically fit well to the data points, but they are 
relatively more common for comer-like sets. Corner-like neighborhoods 
may be better approximated by implicit cubic polynomials, which include 
triplets of planes in their family. Many cases occurring in practice could also 
be improved by approximations within certain sub-families of surfaces, such 
as intersections of planes with circular cylinders or cones. 

Outliers are points that are far from the surface and may be created 
either by mistake or by physical shortcomings of the scanning hardware. 
Outliers may cause trouble in the reconstruction of the surface and are 
preferably removed before they do so. Because the surface is not known, 
outliers can be detected only with indirect methods. A straightforward 
approach for detecting outliers specifies a threshold and removes points 
whose square distances to the planes of best fit exceed that threshold. A 
drawback of this method is that many or all points with edge-like and comer- 
like neighborhoods are likely to be classified as outliers. This method is 
refined by observing that points in such edge or comer regions have near 
points with similarly large square distances to their planes of best fit To 
discriminate between points with and without such near points, the average 
square distance between points and corresponding planes of best fit are 
considered. This average is: 

n i 



-23- 



03/027961 



PCT/USO2/24220 



* The point p, is treated as an outlier rf Its square distance exceeds a constant 
times the average as expressed by: 

(P>,) 2 >C 0 5. 

After computing aii square distances, the outliers are identified by evaluating 
this simple inequality. It can easily be modified by adjusting the constant 
The above-described operations may be treated as preprocessing 
operations that reduce noise and remove outliers from the data points so 
that the quality of surfaces reconstructed therefrom can be improved. 

. Operations to reconstruct three-dimensional (3D) surfaces from 
dense and locally two-dimensionally distributed data points sets, (e.g., point 
clouds) will now be described. As described above, these data point sets _ 
preferably undergo preprocessing operations to reduce noise and remove 
outliers, before reconstruction operations are performed. These 
"reconstruction operations include generating a plurality of stars by locally 
projecting sampled points within a point set onto estimated tangent planes 
and then merging the stars in two-dimensional Delaunay triangulations 
within the estimated tangent planes. The Delaunay triangulations are 
preferably constructed as weighted Delaunay triangulations, however, 
Delaunay triangulations having unweighted vertices may also be 
constructed if less accurate results or less efficient operations are 
acceptable. The reconstruction operations return a two-dimensional 
simplicial complex in which the edges and triangles that share a vertex form 
a portion of an open disk. Such complexes are treated herein as 
pseudomanifolds. 

Although the surface reconstruction operations described herein are 
based on concepts from differential geometry, they may be classified as 
surface meshing operations. The reconstruction operations may be 
arranged in four stages, and reiy on the data point set ScE 3 being densely 
sampled on the surface of an object or otherwise densely generated. These 
four stages include: 
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(A) Finding, for each of a plurality points p, e S, a 
respective subset of near points S, c S, which contains p,; 

(B) Estimating a tangent plane using S h projecting S, 
onto that plane, and constructing the star of p,- in the two- 
dimensional Delaunay triangulation; 

(C) Eliminating edges and triangles that do not belong to 
the stars of all their vertices, and constructing the surface from 
what remains; and 

(D) Performing post-processing by identifying and filling 
holes in the surface created by the elimination of edges and 
triangles during stage C. 

As described above with respect to FIG. 1 , the operations of stage A, 
which include searching for and identifying a respective set of near points S, 
for each of a plurality of points p h can be time-consuming. To improve the 
efficiency of determining a respective set of near points S, for each of a 
plurality of points p h a data structure that exploits the locality of the access is 
used for repeated searching. The near point search operations include . 
storing an input set Sci 3 in an implicit oct-tree and then using a random 
sample R c S to compute a width 2r 0 of cubes having a specified average 
yield m 0 . For each point p, e S, a subset of k Q nearest points that lie in the 
cube of width 2r 0 , centered at p /( is computed. To describe this in 
mathematical terms, the box function Ofr) is defined as the set of points p y e 
S with an /..-distance at most r from point p h where r is a real number. 
Stated alternatively, the box function defines a set of points p y that are 
contained within the closed axis-aligned cube of width 2r centered at p h 
Now using the Euclidean distance, we define N£k) as the set of k points 
closest to p„ including p, itself. The set of near points of p, is then: 

where r 0 is the positive real number computed to achieve a specified 
average yield and kg is a positive integer constant It makes sense to define 
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k 0 a few times the average degree of a vertex in a flat triangulation, which is. 
about six. in a typical case, = 30. 

Operations for constructing an oct-tree will now be described. 
Additional information relating to the construction of oct-trees can be found 
in a two-volume textbook by H. Samat, entitled Spatial Data Structures: 
Quadtrees, Octrees, and Other Heirarchical Methods, Addison-Wesley 
(1989). Upon normalization, the input data point set can be treated as being 
contained in a half-open unit cube, S c [0, 1) 3 . The oct-tree corresponds to 
a recursive decomposition of this half-open cube into eight congruent half- 
open cubes. FIG. 3 illustrates the construction of an oct-tree 30 in two 
dimensions, where a square 32 is recursively decomposed into four smaller 
squares. Each node ft of the oct-tree 30 is assigned an integer address 
that encodes the path from the root to /jj. Three bits per level are required in 
the illustrated example, which allows at most / 0 = 10 levels in order to fit 
each address into a 32 bit word. The oct-tree is constructed from top down, 
with each node being decomposed as long as its level is less than / 0 and the 
number of points in its half-open cube exceeds n Q = 30. 

A leaf is a node that is not decomposed any further. The oct-tree is 
represented by its pre-order sequence of leaves and a parallel sequence of 
points. The ordering provides that the points that lie inside a cube of a node 
are contiguous in the second sequence. Each leaf stores its address and 
the index interval of the points in its cube. Navigation is done through bit 
manipulations of node addresses* The problem of enumerating all points 
inside a cube C will now be considered as an example. Let C(//) be the half- 
open cube of node //, let S(jJ) = Sn C(jJ), and write p for the root of the oct- 
tree. The points in C can be found by calling function POINTS with // = p, 
as illustrated by the following program code: 

void POINTS (C, * * 

if Qjj) q C then return S(jj) endif; 

if C(j/) n C # 0 then 
if y is a leaf then return S(fj) n O endif; 
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forall children v of //do POINTS (C, v) endfor 
endif. ; 

The next operation includes computing a distance r 0 that is small 
enough (but not too small) to estimate differential properties of the surface. 
Two empirically tested constants m 0 = 100 and s 0 - 100 can be used, where 
m 0 represents a desired average yield that may be selected by a user and s 0 
represents the size of the random sample RzS. 

The yield of an axis-aligned cube is defined as the number of points 
of S it contains. For a given real number r > 0, the average yield is 
computed from all axis-aligned cubes of width 2r centered at points in S. 
The real number r 0 is defined as the smallest value of r for which the. 
average yield is at least m^. Computing this value of r is typically time- 
consuming, but it can be estimated quickly by choosing a random sample R 
cSaf size s 0> and computing the average yield m(r) over the cubes 
centered at points in R. The real number r Q is computed as the minimum 
value of r for which m(r) a m a To save time, a small estimatels used as a 
start and then it is improved by first growing and then by shrinking the 
interval, as illustrated by the following program code: 

a- 0.0; b = 0.001; 

while m{b) <m 0 6ob = 2b endwhile; 

for/= 1 to /„ dox = (a + b)/2; 

if m (x) < m 0 then a = x else b = jc endif 

endfor. 

After completing the iteration, r„ = b can be used, which may be a 
little larger than promised, or r 0 can be selected as the Sgm 0 smallest 
distance defined by the points in the sets q(£), over alfp, e R. Important 
operations include (i) computing m(r) and (ii) selecting from a set of • 
distances. The former uses function POINTS illustrated above and the 
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latteris performed using a one-sided version of randomized quicksort, which 
is described more fully hereinbelow. 

Operations to compute the sets S, of near points, for all p,eS will 
now be described. These operations may be the same as those used to 
determine the width of the near point search cube, but speed is more critical 
because the operations to compute sets of near points are applied to more 
points. The value C, is written for the axis-aligned cube of width 2r 0 with 
center p, and p for the root of the oct-tree, as illustrated by the following 
program code: 

fbr/=1to/ido 

N = POINTS(C, p); 

select Sj as subset of k 0 points in N closest to p, 
endfbr. 

A one-sided version of randomized quicksort is used to select points. To 
understand this operation, suppose N is a linear array storing the points in 
an arbitrary order. Quicksort splits the array using a randomly chosen pivot 
point, p. Specifically, function SPLIT rearranges the array so that the points 
to the left of p are closer to p, and the points to the right of p are farther from 
p t The function then returns the position of the pivot Instead of recursively 
sorting on both sides of that position, the recursive sort is performed only on 
one side. Initially, lo and hi are the first and last indices in N and Ac - /c^ 
These operations are illustrated by the following program code: 
integer SELECT (k, lo, hi) 

assert to s lo + 1 * hi; 

if lo = hi then return lo endrfi 

mid = SPLIT (lo t hi); 

if * s mid-lo then SELECT (k, lo, mid- 1) ' 

else rf k = mid - lo + 1 then return mid 

else SELECT(/r - mid + lo - 1 , mid +1 , hi) 
endif. 
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The expected running time for m.input points is O(m), with a reasonably 
small constant of multiplicity. The operation SELECT is called repeatedly, 
with an average size of N being equal to about /fy. Because.the SELECT 
operation is typically executed a large number of times, ft is important that it 
be optimized for this input point size. 

Operations to construct stars, stage B, can be considered as a 
sequence of operations that include: 

(B1) Using points in S,to estimate the normal n, and the 

estimated tangent plane T, :.<x, n,) = 0 at p; 

(B2) Projecting the points in S/to form the set S/ of weighted 

points in T„ which contains points p,'; and 

(B3) Computing the star of p,' in the weighted Delaunay 

triangulation of S,'. 

tt is preferred that operations B1 , B2 and B3 are performed for each point p, 
e S. Each edge of the resulting triangulation can occur in up to two stars 
and each triangle can occur in up to three stars. Moreover, as explained 
more fully hereinbelow, weights are preferably used to counteract distortion 
caused by the projection of the points in S, to form the set S; of weighted 
points in T t 

The operations B1 assume that each set of near points S, is sampled 
from a small and smooth neighborhood of the point pj, and that the points p y 
e Si lie close to the plane tangent to the surface and passing through p t 
Because the surface to be reconstructed is not known, the locations of 
planes tangent to the surface are also not known. Nonetheless, the tangent 
planes may be estimated from respective sets of near points S,. In 
particular, each estimated tangent plane 7", is computed as the two- 
dimensional linear subspace of E 3 that minimizes the sum of square 
distances to the points p y -p„ over all p y e S,. The definitions of T, as the 
estimated tangent plane at p, and n, as the estimated normal at p, are 
illustrated in FIG. 1. 
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The operations B2 include determining weighted points in T h In. 
particular, for each point p, e S„ the projected point pj = (p/, Wj ) is the 
weighted point in 7), where: 

is the orthogonal projection of p-p, onto T t and Wj = -<p y - p /( n,} 2 is the 
negative square distance of that point from T h We use the weight to 
counteract the distortion of the inter-point distances. To explain this, we 
define the weighted square distance of a point x e7", from pj as: 

n/x) = [|x-p;|| 2 -w jf 

which is the Euclidean square distance of x from p t As described in a 
textbook by H. Edelsbrunner, entitled "Geometry and Topology for Mesh 
Generation," Cambridge Univ. Press (2001), the disclosure of which is 
hereby incorporated herein by reference, a weighted Voronoi region of pj is 
the set of points x e7) whose weighted square distance to p/ is no larger 
than any other weighted point, V, = {x el, | t\£k) s n/x), Vj}. It is also the 
intersection of the three-cfimensional (unweighted) Voronoi cell of p, with the 
plane T r The rationale is that as long as T f intersects the Voronoi ceil for p, 
in the same set of facets and edges as the hypothetical surface, the star of 
p! in the weighted Delaunay triangulation (the dual of the weighted Voronoi 
diagram) is the projection of the star of p, in the restricted Delaunay 
triangulation (the dual of the intersection between the three-dimensional 
Voronoi diagram and the hypothetical surface). The concept of a restricted 
Delaunay triangulation is more fully described in the article by H. 
Edelsbrunner and N.R. Shah, entitled Triangulating Topological Spaces," 
Intemat J. Comput. Geom AppL, Vol. 7, pp. 365-378 (1997), the disclosure 
of which is hereby incorporated herein by reference. In other words, the 
stars in the two-dimensional weighted Delaunay triangulations are the 
closest representations to an ideal reconstruction, which is the restricted 
Delaunay triangulation. 
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Understanding the operations B3 for computing the star of p/ in the 
weighted Delaunay triangulation of S/ requires a basic understanding of 
weighted Delaunay triangulations in a piane, which will now be provided. 
The weighted square distance to two weighted points p/ = (p f ", w,) and x' = 
(x,w), can be generalized as: 

n.OO^lp/'-xIf-Wi-w. 

The two weighted points are orthogonal if n/x') = 0, which in geometric 
terms means that the circles centered at p f " and x with radii w, ia and w 1c 
intersect at a right angle. The two weighted points are considered closer 
than orthogonal if n^x') < 0 and farther than orthogonal if ri^x') > 0. In the 
plane, for any three weighted points p,', p/ and p„\ there is a unique 
weighted point x' that is orthogonal to all three. This weighted point x**is the 
orthocenter of p/, p/ and p k \ If SJ is the collection of weighted points in I 2 , 
by duality to the weighted Voronoi diagram, the weighted Delaunay 
triangulation of S; consists of all triangles A"p,"Pk" for which the orthocenter 
of Pi', pf and p k ' is farther than orthogonal to all other weighted points in S}. 
As described at section 1.4 of the aforementioned textbook by H. 
Edelsbruriner, the genericity of the set of weighted points can be simulated 
to avoid ambiguities in the characterization of the weighted Delaunay 
triangulation. An example of a two-dimensional weighted Delaunay 
triangulation 40 is illustrated by FIG. 4, where the shaded star of p/ within 
the weighted Delaunay triangulation of S/ is either an open disk or an open 
half-disk. 

The star of a projected point p/, denoted herein as St p,', consists of 
the point p/ and all edges and triangles in D{ that contain p/ as a vertex. 
The underlying space of the star is the union of the interiors of its simplices, 
as defined by the following relationship: 

!Stp/| = U<KSt P 7intO. 
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Two cases can be distinguished, based on whether the projected point p/ is 
an interior vertex or a boundary vertex. If p/ is an interior vertex of the 
weighted Delaunay triangulation, then the edges and triangles in the star 
alternate and close a ring about the vertex. In this case, the underlying 
space is an open disk. However, if p/ is a boundary vertex of the 
triangulation, then the edges and triangles still alternate, but form only a 
sequence about the vertex. In this case, the underlying space is an open 
half-disk. Both cases are illustrated by the shaded stars in FIG. 4. When 
the projection is not important or reversed, then reference may be made to 
the star of p, instead of the star of p,'. 

The operations B3 include computing the star of p/ in a 
counterclockwise order around p/. As described herein, the distance 
between two weighted points will be computed as the weighted square 
distance between the two weighted points. These operations begin by 
renaming the weighted points such that q Q = p/ ( q 1 is the weighted point 
closest to q 0t and q v q 2t ... f q M is the counterclockwise order of the weighed 
. point around q^ An ambiguity in the order arises when two weighted points 
lie on the same half-line emanating from q 0 and this ambiguity is resolved by 
throwing away the weighted point farther from q^ It is convenient to repeat 
q k = q 1 at the end of the ordering. To distinguish the interior from the 
boundary vertex case, the orientation of point triplets is tested. Using Greek 
letters for the x and y Cartesian coordinates in T h the orientation of the point 
triplet qrs is the sign of the determinant of the matrix I\ where q = (i|r u ijr 2 ) t r 
" (Pu ft) and s = (o^Oi): 

"1 y/ x 
T= 1 p x p 2 

.1 ^1 ^2J 



A positive orientation means that a left-turn is taken at r, coming straight 
from q and going straight to s. The interior vertex case is characterized by 
the property that all triplets q 0 qfl^ have positive orientation. In the 
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10 



15 



boundary vertex case, there is precisely one index j for which does 
not have positive orientation. The vertices are then relabeled so that q >t is 
the first and q, is the last in the counterclockwise ordering around q 0 . In 
contrast to the interior vertex case, the first vertex is not repeated at the end 
of the ordering. 

With these preparations, the operations for constructing the star are 
the same for both interior and boundary vertices and are performed 
incrementally. For each new weighted point, all triangles whose 
orthocenters are closer than orthogonal are removed and then one new 
triangle is added to the star. At any moment, the star is a sequence of 
triangles, as illustrated in FIG. 5. In FIG. 5, if all weights w, are assumed to 
be zero, the new point q s is closer than orthogonal to the orthocenters. of the 
last triangle g 0 g 3 g 4 in the shaded partial star 50. That last triangle is 
removed before q 5 is added as a vertex of the new last triangle qti q 3 q 5 . The 
operations PUSH. POP, TO and isEMPTY are used to manipulate an 
initially empty stack that stores the sequence of vertices on the boundary of 
the partial star 

PUSH 2 (q v <7 2 ); 
for; = 3 to* do 
20 loop r, s- TOP 2 ; 

if TOOCLOSE(g 0 , r, s, qj) then POP 
else PUSH (q ; ); exit 
endif 
forever 

25 endfor. 

The boolean function tooCLOSE tests whether or not the point g y is closer 
than orthogonal to the orthocenters of q* r, and s. Using again Greek 
letters for the coordinates and the weights, the matrix A is defined as- 

30 
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10 



for qf = g 0 , r, s and t = The weighted point t (= (t 1f tj)) is orthogonal to the 
orthocenter of qrs if and only if det A = 0. The triplet qrs has positive 
orientation, by construction, and it follows that the determinant of the upper 
left three-by-three submatrix (in the above matrix A) is positive. Hence, t is 
farther than orthogonal to the orthocenter if and only if det A > 0. The 
following operations are thus defined: 
boolean TOOCLOSE (g, r, s, f) 
return det A < 0. 



The running time of the star construction operations is 0(k log k) for sorting 
plus 0(A) for constructing the star as described above. 

The operations to compute stars may result in stars that share 
triangles and edges with other stars and stars that conflict with other stars, 
15 because the operations are performed independently in the various 

estimated tangent planes. In the event a conflict is present, the edges and 
triangles that are in conflict are eliminated and the remainder of edges and 
triangles that are not in conflict are merged into a single surface description. 

Operations to merge stars, stage C, can be treated as a sequence of 
20 four sub-operations that include: 

(C1) Sorting F, which is a list of all triangles in the stars, and remove 
all triangles that are not in triplicate; 

(02) Connecting the remaining triangles to form a triangulated 
pseudomanifold; 
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(C3) Sorting E, which is a list of edges in stars that do not belong to 
any triangles in the pseudomanifold, and remove ail edges that are not. in 
duplicate; and 

(C4) Adding the remaining edges to the triangulated 
5 pseudomanifold. 

The representation of the pseudomanifold created in operations C2 
and C4 may be referred to as a trist data structure. In the operation C1 to 
sort triangles, each triangle is represented by the ordered triplet of indices of 
10 its vertices, ijk with / < j < k. The sort operation is performed 
lexicographically: 

'i<i\ 

ijk precedes /'/ k'V^i- Vand j<f,or 

i = i',j=f,andk<k\ 

The use of contiguous integers suggests the use of a radix sort operation, 
rather than a comparison-based sorting operation. To describe the radix 
1 5 sort operation, n is defined as the number of points in S and B[1 ..n] is 

defined as a linear array of buckets. Each bucket is initially an empty stack 
of integer pairs. In the first phase, the triangles are spread using the 
smallest of the three vertex indices as the address in B. The operation 
PUSH, Q,k) is used to add the pair jk to the /-th bucket: 
20 forall/do 

forall triangles ppp k e St p, do 
rename indices such that / <j <k, PUSH, (j,k) 
endfor 
endfor. 



25 



A single vertex belongs to at most some constant number of triangles. The 
sorting operation can therefore be finished by calling quicksort for each 
bucket without paying a logarithmic term in comparison to additional phases 
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of a radix sort. After sorting, Bp] contains the triangles ijk lexicographically 
sorted by jk. The triangles that are not in triplicate can thus be eliminated by 
scanning. The operations POP„ TOP,, and ISEMPTY, can be used to 
manipulate the Mh bucket The triangles that occur three times are 
collected in an initially empty list F, which is extended using function ADD. 
The following pseudocode illustrates these operations: 
fbr/=1tofldo 

while not ISEMPTY, do jk = POP;; 
if TOP, =y* then y/r= POP,; 

if TOP, = jk then jk = POP,; ADD (ijk) end if 
endif 
endwhile 
endfor. 

The trist data structure is an array of (unordered) triangles connected 
to each other by adjacency. In the operation C2, each unordered triangle 
.with vertices p„ p p and p k is represented by its six ordered versions. As 
illustrated by FIG. 6, each ordered version of an unordered triangle stores a 
pointer fhext to the next triangle in the ring around the directed edge 60a- 
60f identified by the first two vertices in the ordering. By construction, each 
edge belongs to at most two triangles, which implies that three of the six 
fnext pointers are redundant. We find matching pairs of triangles by sorting 
the 3f (unordered) edges of the ftriangles in F using the radix sort 
operation. Let E 2 be the list of edges obtained after removing the ones that 
occur only once. In this list we represent each edge by the ordered pair of 
vertex indices, ij with / < j t and a tri pointer to the triangle responsible for its 
existence in E 2 . The triangles are finally matched by scanning E 2 . If 2e 2 is 
defined as the length of £ 2( then the following operations can be performed, 
to match up all triangles: 

for / = 1 to ej do 

CONNECT (£ 2 [2 / - 1].tri,E 2 [2 /J.tri) 
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endfor. 

Besides the edges that are shared by two triangles, boundary edges are 
provided in the resulting pseudomanifold. For each boundary edge pp p a 
triangle pp/o is constructed, where w is a new dummy vertex. The dummy 
triangles 70 connecting w to boundary edges are illustrated in FIG. 7. In 
particular, FIG. 7 illustrates five dummy triangles 70a-70e, four of which 
(i.e., 70a-70d) connect u to boundary edges and one of which (i.e., 70e) 
connects o> to a principal edge in the star of p r The boundary edges are 
exactly the ones removed from E 2 . Assuming the boundary edges are 
collected in another list, £,, the dummy triangles are constructed in a single 
scan, as illustrated by the below operations: 

for/sltoe^o 

let/)' be the index pair in E,[/J; 

create pp/o; CONNECT (E,[/].tri, pp/o) 

endfor. 

The operations C3 include adding principal edges. In these 
operations, all. edges that belong to the stars of both their endpoints are 
accepted as members of the pseudomanifold. The edges that belong to 
one or two accepted triangles are already part of the pseudomanifold, but 
the principal edges that belong to no accepted triangles need to be added. 
Each principal edge pp k is represented by the dummy triangle pp^io shown 
as dummy triangle 70e in FIG. 7. All principal edges are found by radix 
sorting the edges in the collection of stars constructed in accordance with 
operations B3. The edges that are not in duplicate are removed along with 
edges that are in E, or in E 2 . The resulting list, £ 0 , stores all principal edges. 
Similar to the boundary edges, the principal edges are' added to the 
pseudomanifold in a single scan, as illustrated by the following operations: 
for/=l toe 0 do 

let /)' be the index pair in £„[/]; create pp/o 
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endfor. . 

It is possible that a problematic case of two 'crossing' principal edges may 
occur. For example, if p h p p p k and p, form a square, then the distortions 
5 caused by the projections of these points into four estimated tangent planes 
may create a contradictory agreement between diagonally opposing 
vertices. This is shown by FIG. 8, where two 'crossing' principal edges 80a 
are shown on the left. As also shown on. the right side of FIG. 8, a similarly 
problematic case may also occur between two accepted triangles 80b. The 

10 case of 'crossing' edges can be avoided relatively inexpensively by 
prohibiting the occurrence of a principal edge pp k for which the two 
neighboring edges pp } and pp t in the star of p, are connected by another 
principal edge, namely pp x . Such cases are again found using a radix sort 
operation and the involved edges are removed from £ 0 before they can be 

1 5 added to the pseudomanifold . 

Operations C4 to establish order and to add the remaining edges to 
the triangulation of the pseudomanifold will now_be described. As described 
above, each star St p, is represented by a list storing the vertices on the 
boundary in order. These lists are used to connect the dummy triangles to 

20 each other. First, boundary and principal edges in the star are identified. 

Each boundary edge either starts or ends a hole when the star is read in a 
counterclockwise order. By definition, a principal edge both ends a hole and 
starts a new one. The buckets B[1 ..n] are used again. The operations 
include storing in Bp] all indices k of vertices that define dummy triangles 

25 PP*to. The buckets are sorted using quicksort, as described above, and 

each index is stored with a pointer to the corresponding dummy triangle. A 
given vertex index j can be located in B[\\ using a binary search operation. If 
j e B[0, then pp } is either a boundary edge that starts a'hole, a boundary 
edge that ends a hole, or a principal edge. The identity of pp y can be 

30 determined by checking the non-null fnext pointers of the dummy triangle 
pp/j>, if any. The boolean operations STARTS and ENDS can be used to 
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express the test The function FIRST, is also used to return the first vertex 
in the Mh star, and the function NEXT, is used to return the next vertex that 
defines a boundary or principal edge. After passing the last such vertex, 
function NEXT, returns null. To reduce the disk and half-disk cases to one, 
the first vertex of a half-disk is repeated at the end of the list, making it 
appear as a disk. These operations are described by the following 
operations: 

for/sltondop^RRST,; 
repeat PjsNEXT; 
If STARTS(py) then assert ENDS(pJ; 

CONNECT (pp/o, ppfii); 
endlf; p y = p, 
until p„ = null 
endfor. 

The most time-consuming part of these operations is the classification of 
edges, which is done using a binary search of the sorted lists SRJ. The 
lengths of these lists is bounded from above by the constant - 1 , which 
implies that searching takes only constant time. The order among the 
dummy triangles is thus established in a time proportional to the number of - 
edges. 

Post-processing operations to identify and fill holes, identified above 
as stage D, include multiple sub-operations that will now be described. 
Holes in a pseudomanifold may be caused by missing data or by conflicts 
between vertex stars. The operations described hereinbelow can be used 
to fill small and simple holes automatically, but large or complicated' holes 
are typically left to a user to manually fill. The operations to fill holes 
include: . . 

(D1 ) Finding holes by identifying boundary cycles; 
(D2) Accepting holes that can be filled automatically; and 
(D3) Filling accepted holes with triangulated polygons. 
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Operations to determine, boundary cycles can be performed efficiently 
because the vertex stars are ordered. However, testing whether or not a 
boundary cycle (e.g., hole) can be filled automatically is typically a more 
complex operation and may depend on the type of shape being 
reconstructed and the application to which the reconstructed surface is to be 
applied. A few conservative heuristics can be used to determine whether or 
not a boundary cycle may be filled automatically. The complexity of the hole 
filling operations is typically proportional to the size of the hole to be filled. 

Operations D1 for identifying boundaries include using the 
classification of edges as principal, boundary, or interior edges, as 
described above. Each boundary edge pp l starts a hole in the star of one 
endpoint and ends a hole in the star of the other endpoint. A directed graph 
H is constructed that represents each principal edge pp y in Kby its two 
directed versions, p^ and pp^ t and each boundary edge that starts a hole in 
St Pi and ends one in St p y as the directed edge pp l . Stated alternatively, 
the directed graph H is a directed version of the link of (o in /C. The set of 
directed edges in H can be partitioned into directed cycles and the boundary 
cycles are identified as respective directed cycles. Let coppj be a triangle in 
the star of o> such that pp l is a directed edge in H. The cycle that contains 
ppj can be traced by following fnext pointers. Function APPEND adds a 
new edge at the end of a list representing the traced cycle, as illustrated by 
the following operations: 

while ppj is unmarked do 

APPEND (pp y ); mark pp y ; a>pp y = cop/vfhext 

endwhile. 

Observe the permutation uppi of copp y used to indicate a different ordered 
version of the triangle. As illustrated in FIG. 9, the cycles are traced such 
that the holes lie locally to the left of the directed edges. We trace cycles 
until all directed edges in H are marked. 
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The operations D2 include accepting holes that can be filled 
automatically, which typically includes only relatively small and simple holes. 
The size of the holes can be measured either geometrically, as the length or 
the diameter, or combinatorially, as the number of directed edges of the 
5 boundary cycle. The size of the holes can be distinguished by introducing a 
size threshold. However, the distinction between simple and complicated 
holes is more intricate, A hole may be treated as a "simple" hole if Its 
boundary cycle can be embedded in a plane. 

Because of the nature of the operations used in the construction of a 
10 cycle, its embedding surrounds an open disk on its left, as shown in FIG. 10. 
This open disk represents the hole to be filled. The indices of the vertices 
that are encountered when a boundary cycle is traversed can be recorded. 
A single vertex may occur an arbitrary number of times, but an alternation of 
the form J.JJ:j„ is not possible. Thus, a proper boundary cycle is a 
1 5 Davenport-Schinzel cycle of order at most two. This can be illustrated 
simply by taking a circle and marking four points with alternating labels, 
i.j. The point / is then connected to point /' by a path outside the circle. 
Similarly, the pointy is connected toy by a path outside the circle. As long 
as the paths are drawn in general position, they will cross an odd number of 
20 times. The circle is then deformed into a boundary cycle so that the ends of 
each path become the same and the paths become closed curves. 
However, two closed curves in general position necessarily cross an even 
number of times, which is a contradiction. 

An operation is used that decides whether or not a cycle Z of indices 
25 contains a forbidden alternation between two indices, /and/ It is 

convenient to cut 2 open to form a sequence. A forbidden alternation 
defines two intervals, one from / to / and the other from j to that overlap but 
neither is contained in the other. Two intervals with this property are 
independent In the first pass, the operation computes all intervals between 
contiguous occurrences of the same index Each interval is represented by 
having its endpoints point to each other. In the second pass, the cycle Z is 



30 
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scanned to see whether there is an independent pair of intervals. If the 
sequence is stored in an array Z[1 ..m], each location / stores the index of a 
vertex, the left endpoint /- of the interval that ends at /, and the right 
endpoint /+ of the interval that starts at /. Both are positions in Z and the 
5 value zero is used to represent the non-existence of such endpoints. The 
operations, which maintain an initially empty stack of currently enclosing 
intervals, are described by the following operations: 
boolean ISSIMPLE (Z) 
for / = 1 to m do 
10 if /-.* 0 then [a, b] = POP; 

if [a, b] # [/-, /] then return FALSE endif 
endif; 

if /+*0 then PUSH ([/, /+]) endif 
endfor, return TRUE. 

15 - 

- Function ISSIMPLE is correct because if all intervals are either disjoint or 
- nested, then the last interval pushed on the stack is the first one to be 
popped. The running time is proportional to the length of the cycle. 

Operations D3 for triangulating holes will now be described. As holes 

20 get larger and more complicated, filling operations that get progressively 
more sophisticated are required. However, if only simple holes are 
considered, the operations may be efficiently performed. These hole filling 
operations create no new vertices and fiH a hole with edges and triangles 
connecting the boundary cycle, Z. By restricting the operations to simple 

25 * holes, these filling edges and triangles form the triangulation of a disk 

whose dual graph is a tree. The triangles that correspond to leaves can be 
referred to as ears. If Z has m > 3 edges, then there are at least two ears. 
The triangulation can thus be built by adding an ear at a time. More 
precisely, the operations find two consecutive edges pp^ and pp k in Z, add 

30 the triangle ppp k and the edge pp k to K, and substitute PAc for p-pj. pp k in Z. 
To produce a reasonable triangulation, the ears are prioritized by the angle 
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at the middle vertex, pj. Because the hole lies locally to the left of pp ]t pp k , 
the angle on that side in the projection is measured onto the estimated 
tangent plane. To avoid duplications, edges pfa that are already part of the 
triangulation surrounding the hole are rejected. The operations use a 
priority queue for the edge pairs, which can be manipulated using the 
functions MINEXTRACT and CHANGEKEY (where p r and pf are used to 
represent the predecessor and successor of p ; in Z): 
while m > 3 do p } = MINEXTRACT; 
rf not MK (pp p^+) then 
add Pj-p^j+ and Pj-Pj+ to K; 
CHANGEKEY 2 (p r , Pj+); /n = /n -1 
endrf 
endwhile; 

Pj = MINEXTRACT; add p f ppf to K. 

The running time for filling a hole bounded by a cycle with m edges is 
0(mlog m). 

The above-described operations for creating and merging stars will 
now be more fully described with respect to FIGS. 1 1-14. In particular, FIG. 
11 illustrates a small set of points, including p 1 and p 2l that define a portion 
of a three-dimensional (3D) surface. In FIG. 12A ( a first set of near points 
S 1 associated with the first point p, and a second set of near points S 2 
associated with the second point p 2 are projected to respective tangent 
planes, shown as T1 and T2. The first set of near points S, and the second 
set of near points S 2 typically contain some points that are common to both 
sets. In the tangent plane T1 , the projection of the first point, p/, represents 
the starting point q 0 for creating the star of p/. Similarly, in the tangent 
plane T2, the projection of the second point, p^, represents the starting 
point q 0 for creating the star of p 2 ': Referring now to FIG. 12B, the first 
triangle in the star of p,' is created by connecting the starting point q 0 to the 
nearest point in the tangent plane T1 , shown as point q 1f by an edge. Then, 



-43- 



03/027961 



PCT/US02/24220 



continuing in a counter-clockwise direction relative to the edge q 0 q 1( the 
starting point q 0 is. connected to the next closest point, shown as q^ The 
points q, and q 2 are also connected by an edge to thereby define the first 
triangle qoq^- The next nearest point in the counter-clockwise direction, q 3l 
is then identified. Referring now to FIG. 12C, an evaluation is then made to 
determine whether the point q 3 is within a circumcircle of the first triangle 
q 0 q 1 q 2 , which is defined as a circle that passes through the vertices of the 
first triangle. If the point q 3 is within a circumcircle of the first triangle qcfliCfc. 
then the point q 3 is added as a vertex of a triangle in the star of p/^qo and 
the point q 2 is discarded. As illustrated by FIG. 12D, the sequence of 
operations illustrated by FIGS. 12B-12C are repeated until the star of p/, 
which frequently has about six (6) triangles, is fully defined as an open disk 
or open haff-disk. 

Referring again to FIG. 1 2B, the first triangle in the star of p 2 ' is 
created by connecting the starting point q 0 to the nearest point in the 
tangent plane T2, shown as point q 1( by an edge. Then, continuing in a 
counter-clockwise direction relative to the edge q 0 q 1( the starting point q 0 is 
connected to the next closest point, shown as q 2 . The points q, and q 2 are 
also connected by an edge to thereby define the first triangle q 0 qiq 2 . The 
next nearest point in the counter-clockwise direction, q 3l Is then identified. 
An evaluation is then made to determine whether the point q 3 in FIG. 12C is 
within a circumcircle of the first triangle q 0 Qiq 2 - |f the P° int <h Is outside a 
circumcircle of the first triangle q^c^-then the point q 3 is added as a vertex 
of a new triangle in the star of p 2 ' and the point q 2 is retained as a vertex. 
As illustrated by FIG. 12D. the sequence of operations illustrated by FIGS. 
12B-12C are repeated until the star of p 2 ' is fully defined. 

As illustrated by FIG. 12D, two triangles A and B, which share a 
common edge, are illustrated with heavy shading. The triangle A in the first 
tangent plane T1 and the triangle A in the second tangent plane T2 have 
vertices that map to the same three points, including the first point p 1( in the 
point set illustrated by FIG. 1 1 . Similarly, the triangle B in the first tangent 
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plane T1 and the triangle B in the second tangent plane f 2 map to the same 
three points, including the first point p,, in the point set illustrated by FIG. 1.1 . 
Accordingly, as illustrated by FIG. 13, the above-described edge and 
triangle sorting operations may be performed to merge the pair of triangles 
A together and merge the pair of triangles B together into a surface 
triangulation, containing the star of p, and the star of p^ 

Preferred operations for creating stars from points projected to a 
plane will now be more fully described with respect to FIGS. 14A-14G. In 
particular, FIGS. 14A-14B illustrate operations that may be used to 
determine how the projected point q 5 is to be connected to other points in 
the star of q 0 as the star, which is shown as an open half-disk, is created. In 
FIGS. 14A-14B, the partial star of q 0 is expanded by evaluating whether the 
next point q 5 in the sequence of points qo.q1.q2.q3.Q4.q5.qo »s inside or outside ■ 
a circumcircle of a triangle defined by the two immediately preceding 
vertices, q 4 and q^ and vertex q^ Because the next point q$ in the 
sequence is outside the illustrated circumcircle, the triangle q 0 q4q 5 , 
illustrated as triangle D, is added and the star of q 0 is defined by triangles A- 
D. However, in FIG. 14C, the next point q 5 in the sequence of points 
q9.q1.q2.q3.q4.q5.qo is inside the circumcircle of q 0 q 3 q 4 . Accordingly, the point 
q 4 is dropped from the sequence and the new triangle q 0 q3q 5 is added as a 
replacement for the triangle q^q^ Thus, as illustrated by FIG. 14D, the 
star of q 0 is defined by triangles A-B and E. 

Referring now to FIG. 14E, the introduction of projection weights 
adds complexity to the operations for determining whether a next projected 
point in a sequence q 0 qiq 2 q3— q^o te to be added as a vertex of a new 
triangle. When projection weights are considered, the preferred operations 
include determining whether or not the next projected point in the sequence, 
illustrated as points "a", tt b" or V, is closer than orthogonal to an orthocenter 
of an orthocircle of triangle q 0 q 3 q4» where each of the vertices q 0) q 3 and q 4 
of the triangle q 0 q 3 q4 has respective non-zero weights associated therewith 
that are function of projection distance (e.g., proportional to (projection 
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-* . - . distance) 2 ). The magnitude of each respective weight is illustrated as 
proportional to a radius of a circle surrounding a respective point As 
illustrated by FIG. 14E, each of the weighted points a a n and V appear to be 
closer than orthogonal (c.to) and point "c" appears to be farther than 
5 orthogonal (f.to). To determine conclusively whether or not a respective 

point is closer than orthogonal, the value of the det A can be determined for 
each weighted point to be evaluated, as described more fully hereinabove. 
For each case where det A < 0, the weighted point t is treated as being 
closer than orthogonal to the orthocenter of the triangle. qrs and for each 

10 case where det A > 0, the weighted point t is treated as . being farther than 
orthogonal. As described more fully in section 1 .4 of the aforementioned 
textbook to H. Edelsbrunner, a geometric technique of symbolically 
perturbing a geometric input can be used advantageously to preclude 
inconsistencies in star building that may occur in cases where a weighted 

15 point is found to be orthogonal to the orthocenter (i.e., det A = 0). 

The operations for determining det A can be understood 
geometrically by considering the two cases illustrated by FIGS. 14F-14G, for 
two points p and q. In both cases, the point q is orthogonal to the 
orthocenter p if: 

20 Hp-qli 2 - -s 2 = 0 

FIG. 14F illustrates the real-real case where ( r 2 , s 2 > 0). In this case, the 
exterior angle a extending between the tangents to the circles surrounding 
points p and q (at the point of intersection between the two circles) can be 
25 evaluated to determine whether or not point q is closer than orthogonal to 
the orthocenter p, as shown by the following relationship: 



<90° f.t.o 
= 90° orthogonal 
> 90° c.t.o 
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FIG. 14G illustrates the real-imaginary case where ( i 2 > 0, s 2 < 0). In this 
case, the interior angle a can be evaluated in accordance with the above 
relationships to determine whether or not point q is closer than orthogonal to 
the orthocenterp. The real-imaginary case presents the possibility of 
having a point q that is physically within an orthocircle centered about the 
orthocenterp, but is nonetheless farther than orthogonal to the orthocenter. 

Referring now to FIG. 15, preferred high-level operations 100 for 
modeling a three-dimensional (3D) surface include determining a plurality of 
stars from a plurality of points in a 3D point set that at least partially 
describes the 3D surface, Block 110, and then merging the plurality of stars 
into a digital model of the 3D surface, Block 120. The operations 100 may 
be performed by a plurality of embodiments. In particular, FIG. 16 illustrates 
operations 200 for modeling a 3D surface that include identifying a subset of 
near points for each of a plurality of points in a 3D point set, Block 210. This 
3D point set may take the form of a point cloud data file, with each data 
point being identified by its Cartesian coordinates. The point cloud data file 
need not possess connectivity information that links or defines relationships 
between respective points therein. The point cloud data files may be 
provided in an ASCII xyz data format by conventional digitizers, including 
those manufactured by Cyberware™, Digibotics™, Laser Design™, 
Steinbichler™, Genex™ and Minolta™, for example. The data files may 
describe 3D surfaces that are closed. As will be understood by those skilled 
in the art of threendimensional geometry, all closed 3D surfaces are either 
star-shaped or non star-shaped. Closed surfaces are "star" shaped if and 
only if there exists at least one point on the interior of the volume bounded 
by the closed surface from which all points on the surface are visible. All 
other surfaces are non star-shaped. Examples of star-shaped surfaces 
include a cube, a sphere and tetrahedron. Examples of non star-shaped 
surfaces include toroids (e.g. donut-shaped) and closed surfaces having 
tunnels and handles. 
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: . * The operations 200 of FIG. 1 6 also include projecting a plurality of 
points from each subset of near points to a respective estimated tangent 
plane, Block 220, and then determining a respective star from each set of 
projected points on an estimated tangent plane, Block 230. As described 
above, the operations to determine a star from a set of projected points 
preferably includes assigning weights based on projection distance to at 
least a plurality of the projected points and using these weights to determine 
whether a projected point is to be a vertex of a triangle within the star or is 
to be discarded. The stars associated with each of a plurality of tangent 
planes are then merged into a 3D surface triangulation of the 3D point set, 
Block 240. As illustrated by FIG. 17, these operations to merge stars, Block 
240, preferably include sorting triangles within the plurality of stars and . 
removing those sorted triangles that are not in triplicate, Block 242. The 
remaining triangles that are not removed are then connected to define a . 
triangulated pseudomanifold, Block 244. Operations may also be performed 
to sort those edges within the plurality of stars that do not belong to any of 
the triangles in the triangulated pseudomanifold, and remove those sorted 
edges that are not in duplicate, Block 246. The remaining edges that have . 
not been removed are then added to the triangulated pseudomanifold, Block 
248. 

Preferred operations 300 to model 3D surfaces may also include 
preprocessing operations that improve the quality of the 3D point set by, 
among other things, reducing noise and removing outliers, as illustrated by 
FIG. 18. The preprocessing operations may include determining a 
respective set of near points for each of a plurality of points in a 3D point 
set, Block 310, and then fitting each set of near points with a respective 
approximating surface, Block 320. A denoising operation may then be 
performed by moving each of the plurality of points in the 3D point set to the 
approximating surface that is associated with its respective set of near 
points, Block 330. 
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More detailed embodiments of these preprocessing operations are 
Hiifetrated by FIG. 19. In particular, operations 400 to model 3D surfaces 
may include determining an estimated normal for each of a plurality of 
points in a 3D point set that at least partially describes the 3D surface, Block 
410, and then estimating principal curvature directions on the 3D surface by 
evaluating a differential structure of the estimated normals associated with 
the plurality of points, Block 420. A local neighborhood of each of the 
plurality of points is also classified in terms of its shape characteristic, Block 
430. Based on its classification, an approximating surface is determined for 
each of the local neighborhoods, Block 440. The quality of the 3D point set 
is then improved by moving each of the plurality of points to a respective 
approximating surface that is associated with a local neighborhood of the 
respective point, Block 450. 

Accordingly, the operations to model 3D surfaces may include 
operations to reduce noise within a 3D point set and then use the improved 
3D point set to create and merge stars into a triangulated model of the 3D 
surface. Thus, as illustrated by FIG. 20, these modeling operations 500 
may include determining a respective set of near points for each of a first 
plurality of points in a first 3D point set, Block 510, and then fitting each set 
of near points with a respective approximating surface, Block 520. An 
operation is then performed to generate a second 3D point set by denoising 
the first 3D point set, Block 530. Operations are then performed to 
determine a plurality stars from a second plurality of points in the second 3D 
point set, Block 540. The stars are then merged into a digital model of the 
3D surface using preferred sorting and removal operations, Block 550. 

Referring now to FIG. 21, a general hardware description of a custom 
CAD/CAM workstation 600 is illustrated as comprising, among other things, 
software and hardware components that perform the operations described 
above, including those illustrated by FIGS. 15-20. The workstation 600 
preferably includes a computer-aided design tool 620 that may accept a 
point cloud data set via a file 604, a scanner 606, data bus 602 or other 
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conventional means. A display 610 and a three-dimensional printer 608 are 
also preferably provided to assist in performing the operations of the present 
invention. The hardware design of the above described components 604, 
606, 608 and 610 is well known to those having skill in the art and need not 
be described further herein. The workstation 600 preferably comprises a 
computer-readable storage medium having computer-readable program 
code embodied in the medium. This computer-readable program code is 
readable by one or more processors within the workstation 600 and tangibly 
embodies a program of instructions executable by the processor to perform 
the operations described herein. 

In the drawings and specification, there have been disclosed typical 
preferred embodiments of the invention and, although specific tertns are 
employed, they are used in a generic and descriptive sense only and not for 
purposes of limitation, the scope of the invention being set forth in the 
following claims. These claims include method, apparatus and computer 
program product claims. The method claims include recitations that may 
also be provided as. recitations within apparatus and computer program 
product claims. In particular, the method claims may recite steps that can 
be treated as operations performed by apparatus and/or instructions and 
program code associated with computer program products. 
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THAT WHICH IS CLAIMED IS: 

1. A method of modeling a three-dimensional (3D) surface, comprising 
the steps of: 

determining a plurality of stars from a plurality of points ft in a 3D point 
set S that at least partially describes the 3D surface, by projecting the 
plurality of points p, onto planes T t that are each estimated to be tangent 
about a respective one of the plurality of points p * and 

merging the plurality of stars into a digital model of the 3D surface. 

2. The method of Claim 1, wherein said determining step comprises 
identifying a respective subset of near points S, for each of the plurality of 
points p,. 

3. The method of Claim 2, wherein said determining step comprises 
projecting a plurality of points p l in each subset of near points S; to a 

" respective estimated tangent plane T,. 

4. The method of Claim 3 ( wherein said determining step comprises 
determining for each of a plurality of estimated tangent planes, T„ a star of 
the projection of a respective point p, onto the estimated tangent plane T,. 

5. The method of Claim 4, wherein the star of the projection of a 
respective point p, onto the estimated tangent plane T, constitutes a two- 
dimensional (2D) Delaunay trianguiation. 

6. The method of Claim 4, wherein the star of the projection of a 
respective point p, onto the estimated tangent plane Tj constitutes a two- 
dimensional (2D) weighted Delaunay trianguiation. 
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7. The method of Claim 6, wherein a distance between at least one 
point p, and its orthogonal projection onto a respective estimated tangent 
plane T, is non-zero. 

8. The method of Claim 1 , wherein said determining step comprises 
identifying a respective subset of near points S, for each of the plurality of 
points Pj by storing the 3D point set S in an pet-tree. 

9. The method of Claim 2, where: 

and CDj(r 0 ) is defined as the set of points Pj e S with an /..-distance at most r 0 
from Pj and N^) is the set of ko points that are closest in Euclidean 
5 distance to p-„ including p ; itself, and is a positive integer. 

10. The method of Claim 1, wherein said determining step comprises 
identifying a respective subset of near points S, for each of the plurality of 
points ^ by determining a width of a near point search cube using a random 
sample R c S and identifying those points in S that are within a respective 

5 near point search cube that is centered about each of the plurality of points 
Pi- 
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1 1 . The method of Claim 1 , wherein said determining step comprises 
identifying a respective subset of near points S, for each of the plurality of 
points Pi by: - 

storing the 3D point set S in an oct-tree; 

determining a width 2r 0 of a near point search cube using a random 
sample RcS, where r 0 is a positive real number; and then, for each of the 
plurality of points p„ 

determining a subset of kj, points that are closest in Euclidean distance 
to ft and selecting from the subset all points that are also within a respective 
near point search cube that is centered about a corresponding point p, and 
has a width equal to 2r 0( where k<j is a positive integer. 

12. The method of Claim 1, wherein said determining step comprises 
identifying a respective subset of near points Si for each of the plurality of 
points p, by: 

storing the 3D point set S in an oct-tree; 

determining a width 2r 0 of a near point search cube using a random 
sample RcS, where r 0 is a positive real number that equals a minimum 
value of r for which an average of a yield is greater than or equal to m 0( 
where m 0 is a positive constant and the yield equals the number of points in 
S that are within a near point search cube of width 2r centered about a 
respective point in the random sample R; and then, for each of the plurality 
of points ft in the 3D point set S, 

determining a subset of ko points that are closest in Euclidean distance 
to ft and selecting from the subset all points that are also within a respective 
near point search cube that is centered about a corresponding point p and 
has a width equal to 2r 0 , where ko is a positive integer. 

1 3. The method of Claim 1 1 , where m 0 is about equal to 100 and ko is 
about equal to 30. 



-53- 



03/027961 



PCT/US02/24220 



14. A method of modeling a three-dimensional. (3D) surface, comprising 
the steps of: 

determining a plurality of stars from a plurality of points ft in a 3D point 
set S that at least partially describes the 3D surface, by projecting each of 
the plurality of points p, onto a respective plane; and 

merging the plurality of stars into a model of the 3D surface by 
eliminating edges and triangles from the plurality of stars that are in conflict 
and merging nonconflicting edges and triangles from the plurality of stars 
into a 3D surface triangulation. 

15. The method of Claim 14, wherein said merging step comprises: 
sorting triangles within the plurality of stars and removing those sorted 

triangles that are not in triplicate; 

connecting the sorted triangles that have not been removed to define a 
triangulated pseudomanifold as a two-dimensional simplicial complex in 
which edges and triangles of a star that share a vertex form a portion of an 
open disk; 

sorting edges within the plurality of stars that do not belong to any of the 
triangles in the triangulated pseudomanifold and removing those sorted 
edges that are not in duplicate; and 

adding the sorted edges that have not been removed to the triangulated 
pseudomanifold. 

16. The method of Claim 14, wherein said merging step comprises: 
sorting triangles within the plurality of stars and removing those sorted 

triangles that are not in triplicate. 

17. The method of Claim 16, wherein said merging step further 
comprises: 

connecting the sorted triangles that have not been removed to define a 
triangulated pseudomanifold as a two-dimensional simplicial complex. 
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18. The method of Claim 17, wherein said merging step further 
comprises: 

sorting edges within the plurality of stars that do not belong to any of the 
triangles in the triangulated pseudomanifold and removing those sorted 
5 edges that are not in duplicate; and 

adding the sorted edges that have not been removed to the triangulated 
pseudomanifold. 

19. A method of modeling a three-dimensional (3D) surface, comprising 
the steps of: 

determining a plurality of triangulated neighborhoods from a plurality of 
points in a 3D point set that at least partially describes the 3D surface, by 
5 projecting each of the plurality of points and one or more neighboring points 
in the 3D point set to a respective plane; and 

merging the plurality of triangulated neighborhoods into a digital model 
of the 3D surface. 

20. The method of Claim 19, wherein said determining step comprises 
projecting each of the plurality of points and one or more neighboring points 
in the 3D point set to a respective tangent plane. 

21. A method of modeling a surface of an object, comprising the steps 

of: 

projecting each point and corresponding set of one or more neighboring 
points in a point set to a respective plane; 
5 determining a star for each plane; and 

merging the stars into a surface triangulation. 
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22. The method of Claim 21 , wherein said projecting step comprises 
projecting each point and corresponding set of one or more neighboring 
points in a point set to a respective plane that is estimated to be tangent 
about the point 

. 23. The method of Claim 21 , wherein said step of determining a star for 
each plane comprises determining a star having vertices that are weighted 
as a function of projection distance for each plane. 

24. A method of modeling a three-dimensional (3D) surface, comprising 
the steps of: 

determining an estimated normal for each of a plurality of points in a 3D 
point set that at least partially describes the 3D surface; 
5 evaluating a differential structure of the estimated normals associated 

with the plurality of points to estimate principal curvature directions on the 
3D surface and classify a respective local neighborhood of each of the 
plurality of points in terms of its shape characteristic; 

determining a respective approximating surface for each of the local 
10 neighborhoods; and 

denoising the 3D point set by moving each of the plurality of points to a 
respective approximating surface that is associated with a local 
neighborhood of the respective point. 
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25. A method of modeling a three-dimensional (3D) surface, comprising 
the steps of: 

determining a respective set of near points S, for each of a plurality of 
points p, in a 3D point set S that at least partially describes the 3D surface, 
where SjcS; 

determining a normal bundle for the 3D point set S by determining a 
respective plane h, of best fit for each of the sets of near points S, and a 
respective normal n, for each of the planes h, of best fit; and 

determining a respective approximating surface for each of the sets of 
near points Sj using the normal bundle to estimate respective principal 
curvature directions for each of the sets of near points S,. 

26. A method of modeling a three-dimensional (3D) surface, comprising 
the steps of: 

determining a respective set of near points S, for each of a first plurality 
of points p1 ( in a first 3D point set S1 that at least partially describes the 3D 
surface, where S,c S1; 

fitting each set of near points S, with a respective approximating surface; 

denoising the first 3D point set S1 into a second 3D point set S2 by 
moving at least some of the first plurality of points p1, in the first 3D point set 
SI to the approximating surfaces associated with their respective sets of 
near points S,; 

determining a plurality of stars from a second plurality of points p2 s in the 
second 3D point set S2, by projecting the second plurality of points p2, onto 
planes T, that are estimated to be tangent about a respective one of the 
second plurality of points p2,; and 

merging the plurality of stars into a digital model of the 3D surface. 
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,; ( . 27. The method of Claim 26, wherein said step of fitting each set of 

near points S; with a respective approximating surface comprises fitting a 
first set of near points S, with a first approximating surface by: 

determining respective planes hj of best fit for each of a plurality of 
5 points pj in the first set of near points and 

determining an estimated normal rij for each of the points Pj as a normal 
of its respective plane h, of best fit. 

28. The method of Claim 26, wherein said step of fitting each set of 
near points S, with a respective approximating surface comprises fitting a 
first set of near points with a first approximating surface by: 

determining respective planes hj of best fit for each of a plurality of 
5 points ^ in the first set of near points S,; 

determining an estimated normal n } for each of the points ^ as a normal 
of its respective plane h| of best fit; and 

classifying the first set of near points S, in terms of its shape 
characteristic, by determining estimates of principal curvature directions for 
10 a point p1, from a plurality of the estimated normals 

29. The method of Claim 26, wherein said step of fitting each set of 
near points S t with a respective approximating surface comprises classifying 
a shape characteristic of each set of near points Sj as plane-like and/or * 
edge-like and/or corner-like. 
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30. A method of modeling a three-dimensional (3D) surface, comprising 
the steps of: 

determining a respective set of near points for each of a plurality of 
points in a 3D point set that at least partially describes the 3D surface; 

fitting each set of near points with a respective approximating surface 
that is a selected from a group consisting of cylinders and quadratic and/or 
cubic surfaces; and 

denoising the 3D point set by moving each of the plurality of points in 
the 3D point set to the approximating surface associated with its respective 
set of near points. 

31. A method of modeling a three-dimensional (3D) surface, comprising 
the steps of: 

determining a respective set of near points for each of a first plurality of 
points in a 3D point set that at least partially describes the 3D surface; and 
determining an estimated normal for each of the first plurality of points 

by: 

determining a respective plane of best fit for each of the sets of 
near points; and 

determining a normal for each of the planes of best fit. 

, 32. A method of modeling a three-dimensional (3D) surface, comprising 
the steps of: 

determining a respective set of near points for each of a first plurality of 
points in a 3D point set that at least partially describes the 3D surface; 

determining a normal bundle by determining a respective plane of best 
fit for each of the sets of near points and a normal for each of the planes of 
best fit; and 

determining from the normal bundle at least one respective principal 
curvature direction for each of the sets of near points. 
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33. A method of modeling a three-dimensional (3D) surface, comprising 
the step of: 

denoising a 3D point set that at least partially describes the 3D surface 

by: 

classifying a first neighborhood of points in the 3D point set S1 
using a mass distribution matrix of the first neighborhood of points 
to estimate first normals associated with the first neighborhood of 
points and a normal distribution matrix of the first normals to 
estimate principal curvature directions; 

fitting an approximate surface, which is selected from a group 
consisting of cylindrical, quadratic and cubic surfaces-, to the first 
neighborhood of points; and 

moving at least one point in the first neighborhood of points to 
the approximate surface. 

34. A method of modeling a three-dimensional (3D) surface, comprising 
the steps of: 

identifying a respective subset of near points for each of a plurality of 
points in a 3D point set that at least partially describes the 3D surface by: 

determining dimensions of a near point search space using a 
random sample of the 3D point set; and 

selecting, for each" of the plurality of points, a respective set of 
points in the 3D point set that are within a respective near point 
search space that is oriented about a respective one of the plurality 
of points; 

determining a plurality of stars from the plurality of points in the 3D point 
set by projecting the points in each subset of near points to a respective 
plane; and 

merging the plurality of stars into a digital model of the 3D surface. 
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35. The method of Claim 34, wherein each near point search space 
comprises a space selected from the group consisting of cubes, elipsoids, 
spheres and parallelpipeds. 

36. The method of Claim 34, wherein said step of determining a 
plurality of stars comprises determining a plurality of stars from the plurality 
of points in the 3D point set by projecting the points in each subset of near 
points to a respective estimated tangent plane. 

37. A method of reconstructing a surface of an object from a three- 
dimensional (3D) point cloud data set S derived from scanning the object, 
comprising the steps of: 

determining a respective subset of near points S, c S for each of a 
plurality of points p ( e S; 

estimating a tangent plane T, for each subset of near points S t ; 

projecting each subset of near points S t onto its respective tangent 
plane T,; 

constructing a respective star of each of the plurality of points p, from 
the projected points on each of the tangent planes T,; 

merging the stars associated with the tangent planes T, into a 3D model 
of the surface by eliminating edges and triangles from the stars that are in 
conflict and merging nonconfiicting edges and triangles from the stars into a 
3D surface triangulation; and 

filling one or more holes in the 3D surface triangulation. 
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38. A method of reconstructing a surface of an object from a three- 
dimensional (3D) point cloud, comprising the steps of: 

determining for each of a first plurality of points in the point cloud a 
respective approximating surface that fits the point's neighborhood; 
5 moving each of the first plurality of points to its respective approximating 

surface; 

determining an estimated tangent plane for each of a second plurality of 
points that have been moved to a respective approximating surface; 

projecting each of the second plurality of points and points in their 
10 respective neighborhoods to a respective one of the estimated tangent 
planes; 

constructing 'stars from points projected to the estimated tangent planes; 

and 

merging the stars into a surface triangulation. 

39. The method of Claim 38, further comprising the step of filling one or 
more holes in the surface triangulation. 

40. A method of reconstructing a surface of an object from a three- 
dimensional (3D) point cloud, comprising the steps of. 

denoising the point cloud; 

determining an estimated tangent plane for each of a plurality of points 
5 in the denoised point cloud; 

projecting each of the plurality of points and other points in its respective 
neighborhood to a respective one of the estimated tangent planes; 

constructing stars from points projected to the estimated tangent planes; 

and 

10 merging the stars into a surface triangulation. 
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41 . The method of Claim 40, further comprising the step of filling one or 
more holes in the surface triangulation by: 

constructing a directed graph that represents each principal edge of the 
surface triangulation by its two directed versions and each boundary edge 
5 ■ as a single directed edge; and 

identifying a boundary cycle of at least one hole by partitioning the 
directed graph into directed cycles. 

42. The method of Claim 41, wherein said step of identifying a 
boundary cycle is followed, by the step of identifying simple holes. 

43. The method of Claim 42, wherein said step of identifying simple 
holes comprises identifying holes having index cycles that are Dayenport- 
Schinzel cycles of order less than three. 

44. .A method of denoising a three-dimensional (3D) point set, 
comprising the steps of: 

estimating directions of a local collection of normals associated with a 
local collection of data points in the 3D point set by determining 
eigenvectors of a mass distribution matrix of the local collection of data 
points; and 

estimating directions of curvature associated with the local collection of 
data points by determining eigenvectors of a normal distribution matrix of 
the local collection of normals. 
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;u 45. A method of modeling a three-dimensional (3D) surface, comprising 

the steps of: 

moving each of a plurality of first points in a point set that at least 
partially describes the 3D surface to a respective approximating surface that 
5 is derived by evaluating a respective first point and a plurality of its 
neighboring points in the point set; 

projecting at least one of the first points, which has been moved to an 
approximating surface, to a first plane that is estimated to be tangent about 
the at least one of the first points; 
10 projecting a plurality of points in a neighborhood of the at least one of 

the first points to the first plane; and 

generating a star from a plurality of projected points on the first plane. 

46. A method of modeling a three-dimensional (3D) surface, comprising 
. the step of: 

determining a star of a first point in a 3D point set that at least partially, 
describes the 3D surface, by: 
5 . . projecting the first point and second, third and fourth points in a 
neighborhood of the first point to a plane; 

assigning respective weights to each of the second, third and 
fourth points that are based on projection distance; and 

evaluating whether the projection of the fourth point is within an 
1 0 orthocircle defined by a triangle having projections of the first, 

second and third points as vertices. 
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47. A method of modeling a three-dimensional (3D) surface, comprising 
the step of: 

determining a first star of a first point in a 3D point set that at least 
partially describes the 3D surface, by: 

projecting the first point and first near points in a neighborhood 
of the first point to a first plane that is estimated to be tangent to the 
first point 

assigning respective weights to. each of the projected first near 
points that are based on distances between the first near points and 
the projected first near points; and 

connecting the projected first point and at least some of the 
projected first near points with triangles that share the projected first 
point as a vertex, by evaluating whether a next projected near point 
in a first sequence of projected first near points is closer than 
orthogonal to an orthocenter of a triangle having the projected first 
point and two of the projected first near points as vertices. 
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48. The method of Claim 47, further comprising the step of: 
determining a second star of a second point in the 3D point set by: 

projecting the second point and second near points in a 
neighborhood of the second point to a second plane that is 
5 estimated to be tangent to the second point 

assigning respective weights to each of the projected second 
near points that are based on distances between the second near 
points and the projected second near points; and 

connecting the projected second point and at least some of the 
10 projected second near points with triangles that share the projected 

second point as a vertex, by evaluating whether a next projected 
near point in a second sequence of projected second near points is 
closer than orthogonal to an orthocenter of a triangle having the 
projected second point and two of the. projected second near points 
15 as vertices. 

49. A method of modeling a three-dimensional (3D) surface, comprising 
the step of: 

sequentially connecting a neighborhood of projected points on a plane 
to a first projected point on the plane by evaluating whether at least one 
5 projected point in the neighborhood of projected points is closer than 
orthogonal to an orthocircle defined by a triangle containing the first 
projected point and two projected points in the neighborhood of projected 
points as vertices, with the neighborhood of projected points having weights 
associated therewith that are each a function of a projection distance 
10 between a respective projected point in the neighborhood of projected 
points and a corresponding point in a 3D point set that at least partially 
describes the 3D surface. 

50. The method of Claim 49, wherein the connected neighborhood of 
projected points constitutes a weighted Delaunay triangulation. 
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51 . The method of Claim 49, wherein the plane is estimated to be 
tangent to the 3D surface. 

52. A method of modeling a three-dimensional (3D) surface, comprising 
the steps of: 

projecting a first point in a 3D point set that at least partially describes 
the surface and a set of points in a neighborhood of the first point to a plane 
that is estimated to be tangent to the surface at the first point; and 

creating a weighted Delaunay trianguiation comprising triangles that 
share a projection of the first point in the plane as a vertex and include at 
least some of the projections of the set of points in the neighborhood of the 
first point as vertices that are weighted as a function of projection distance. 

53. The method of Claim 52, wherein the vertices are weighted as a 
function of projection distance squared. 

54. An apparatus for modeling a three-dimensional (3D) surface, 
comprising: 

means for determining a plurality of stars from a plurality of points Pj in a 
3D point set S that at least partially describes the 3D surface, by projecting 
the plurality of points p, onto planes Tj that are each estimated to be tangent 
about a respective one of the plurality of points p,; and :- 

means for merging the plurality of stars into a digital model of the 3D 
surface. 

55. The apparatus of Claim 54, wherein said determining means 
comprises means for identifying a respective subset of near points S, for 
each of the plurality of points p,. 
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" 56. The apparatus of Claim 55, wherein said determining means 
comprises means for projecting a plurality of points p, in each subset of near 
points S, to a respective estimated tangent plane T,. 

57. The apparatus of Claim 56, wherein said determining means 
comprises means for determining for each of a plurality of estimated tangent 
planes, T„ a star of the projection of a respective point p, onto the estimated 
tangent plane T,. 

58. The apparatus of Claim 57, wherein the star of the projection of a 
respective point p, onto the estimated tangent plane T| constitutes a two- 
dimensional (2D) weighted Delaunay triangulatiori. 

59. An apparatus for modeling a three-dimensional (3D) surface, 
comprising: 

means for determining a plurality of stars from a plurality of points Pj in a 
3D point set S that at least partially describes the 3D surface, by projecting 
5 each of the plurality of points p, onto a respective plane; and 

means for merging the plurality of stars into a model of the 3D surface 
by eliminating edges and triangles from the plurality of stars that are in 
conflict and merging nonconflicting edges and triangles from the plurality of 
stars into a 3D surface triangulation. 

60. The apparatus of Claim 59, wherein said merging means comprises 
means for sorting triangles within the plurality of stars and removing those 
sorted triangles that are not in triplicate. 
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61. The apparatus of Claim 60, wherein said merging means further 
comprises: 

means for connecting the sorted triangles that have not been removed 
to define a triangulated pseudomanifold as a two-dimensional simplicial 
5 complex. 

62. The apparatus of Claim 61, wherein said merging means further 
comprises: 

means for sorting edges within the plurality of stars that do not belong to 
any of the triangles in the triangulated pseudomanifold and removing those 
5 sorted edges that are not in duplicate; and 

means for adding the sorted edges that have not been removed to the 
triangulated pseudomanifold. 

63. A computer program product that models three-dimensional (3D) 
surfaces and comprises a computer-readable storage medium having 
computer-readable program code embodied in said medium, said computer- 

• readable program code comprising: 
5 computer-readable program code that determines a plurality of stars 

from a plurality of points p, in a 3D point set S that at least partially describes 
the 3D surface, by projecting each of the plurality of points p t onto a 
respective plane; and 

computer-readable program code that merges the plurality of stars into 
0 a model of the 3D surface by eliminating edges and triangles from the 
plurality of stars that are in conflict and merging nonconflicting edges and 
triangles from the plurality of stars into a 3D surface triangulation. 



-69- 



03/027961 



PCTYUS02/24220 



64. The computer program product of Claim 63, wherein said 
computer-readable program code that merges the plurality of stars 
comprises: 

computer-readable program code that sorts triangles within the plurality 
of stars-and removes those sorted triangles that are not in triplicate; 

computer-readable program code that connects the sorted triangles that 
have not been removed to define a triangulated pseudomanifold as a two- 
dimensional simplicial complex in which edges and triangles of a star that 
share a vertex form a portion of an open disk; 

computer-readable program code that sorts edges within the plurality of 
stars that do not belong to any of the triangles in the trianigulated . 
pseudomanifold and removes those sorted edges that are not in duplicate; 
and 

computer-readable program code that adds the sorted edges that have 
not been removed to the triangulated pseudomanifold. 

65. A computer program product that models three-dimensional (3D) 
surfaces and comprises a computer-readable storage medium having 
computer-readable program code embodied in said medium, said computer- 
readable program code comprising: 

computer-readable program code that projects a first point in a 3D point 
set that at least partially describes a 3D surface and a set of points in a 
neighborhood of the first point to a plane that is estimated to be tangent to 
the 3D surface at the first point; and 

computer-readable program code that creates a weighted Delaunay 
triangulation comprising triangles that share a projection of the first point in 
the plane as a vertex and include at least some of the projections of the set 
of points in the neighborhood of the first point as vertices that are weighted 
as a function of projection distance. 
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66. A method of modeling a three-dimensional (3D) surface, comprising 
the steps of: 

projecting a first point in a 3D point set that at least partially describes 
the surface and a set of points in a neighborhood of the first point to a plane 
that is estimated to be tangent to the surface at the first point; and 

creating a weighted Delaunay triangulation comprising triangles that 
share a projection of the first point in the plane as a vertex and include at . 
least some of the projections of the set of points in the neighborhood of the 
first point as vertices that are weighted as a function of projection distance, 
by evaluating whether or not one or more of the projections of the set of 
points in the neighborhood of the first point are closer than orthogonal to an 
orthocenter of a first triangle in the weighted Delaunay triangulation. 

67. The method of Claim 66, wherein each of at least a plurality the 
vertices is weighted as a function of projection distance squared. 

68. The method of Claim 66, wherein said creating step comprises 
evaluating a matrix containing coordinates of the vertices of the first triangle 
as entries therein. 

69. The method of Claim 68, wherein said creating step comprises 
computing a determinant of the matrix. 

70. The method of Claim 69, wherein the matrix is a 4x4 matrix; and 
wherein at least some of the entries in the matrix are functionally dependent 
on the weights associated with the vertices of the first triangle. 
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71. A surface modeling apparatus, comprising: 

means for performing the method of any one of claims 1-53 and 66-70. 

72. A computer program product readable by a machine and tangibly 
embodying a program of instructions executable by the machine to perform 
the method of any one of claims 1-53 and 66-70. 
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